Main conclusions: Statistical power higher for emergence rate than for number of adults Emergence rate is density dependent with optimum around 100eggs per tube Variation in the number of eggs between 50 and 100 has more effect in emeergence rates than variation between 150 and 250 The SA model should include the number of eggs as a covariate or a factor wit an egg score to increase power

1 Load packages

library(lattice)
library(lme4)
## Loading required package: Matrix

2 Performance

2.1 Load data

data <- read.table(file=here::here("data", "DATACOMPLET_PERF.csv"), header=TRUE, sep=",")

head(data)
##   Generation Experiment Original_environment   Population     Date_P   Date_C_O
## 1         G0 Plasticity                  WT3          WT3 18/09/2018 19/09/2018
## 2         G0 Plasticity                  WT3          WT3 18/09/2018 19/09/2018
## 3         G0 Plasticity                  WT3          WT3 18/09/2018 19/09/2018
## 4         G0 Plasticity                  WT3          WT3 18/09/2018 19/09/2018
## 5         G0 Plasticity                  WT3          WT3 18/09/2018 19/09/2018
## 6         G0 Plasticity           Blackberry Blackberry45    3/10/18    4/10/18
##     Date_C_A Row Column Rack Test_environment Nb_eggs Obs_O Nb_adults Obs_A
## 1    4/10/18  L1     C1    1       Strawberry      41    LO         2    LO
## 2    4/10/18  L1     C2    1       Blackberry      85    LO        32    LO
## 3    4/10/18  L1     C3    1               GF      63    LO        14    LO
## 4    4/10/18  L1     C4    1           Cherry      77    LO        28    LO
## 5    4/10/18  L1     C5    1            Grape      38    LO        19    LO
## 6 19/10/2018  L1     C1    1       Blackberry       0    LO         0    LO
str(data)
## 'data.frame':    2150 obs. of  15 variables:
##  $ Generation          : Factor w/ 2 levels "G0","G2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Experiment          : Factor w/ 2 levels "Adaptation","Plasticity": 2 2 2 2 2 2 2 2 2 2 ...
##  $ Original_environment: Factor w/ 4 levels "Blackberry","Cherry",..: 4 4 4 4 4 1 1 1 1 1 ...
##  $ Population          : Factor w/ 26 levels "Blackberry31",..: 26 26 26 26 26 13 13 13 13 13 ...
##  $ Date_P              : Factor w/ 33 levels "1/10/18","10/7/18",..: 11 11 11 11 11 26 26 26 26 26 ...
##  $ Date_C_O            : Factor w/ 33 levels "1/8/18","11/7/18",..: 11 11 11 11 11 27 27 27 27 27 ...
##  $ Date_C_A            : Factor w/ 33 levels "1/12/18","11/10/18",..: 28 28 28 28 28 16 16 16 16 16 ...
##  $ Row                 : Factor w/ 10 levels "L1","L10","L2",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Column              : Factor w/ 10 levels "C1","C10","C2",..: 1 3 4 5 6 1 3 4 5 6 ...
##  $ Rack                : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Test_environment    : Factor w/ 5 levels "Blackberry","Cherry",..: 5 1 3 2 4 1 3 2 4 5 ...
##  $ Nb_eggs             : int  41 85 63 77 38 0 0 3 1 0 ...
##  $ Obs_O               : Factor w/ 6 levels "","AA","LO","NR",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ Nb_adults           : int  2 32 14 28 19 0 0 2 2 0 ...
##  $ Obs_A               : Factor w/ 6 levels "","AA","CD","JL",..: 5 5 5 5 5 5 5 5 5 5 ...
data$Generation
##    [1] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##   [25] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##   [49] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##   [73] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##   [97] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [121] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [145] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [169] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [193] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [217] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [241] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [265] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [289] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [313] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [337] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [361] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [385] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [409] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [433] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [457] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [481] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [505] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [529] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [553] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [577] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [601] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [625] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [649] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [673] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [697] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [721] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [745] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [769] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [793] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [817] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [841] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [865] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [889] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [913] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [937] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [961] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
##  [985] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [1009] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [1033] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [1057] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G2 G2 G2 G2 G2 G2 G2
## [1081] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1105] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1129] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1153] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1177] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1201] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1225] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1249] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1273] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1297] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1321] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1345] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1369] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1393] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1417] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1441] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1465] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1489] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1513] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1537] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1561] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1585] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1609] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1633] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1657] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1681] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1705] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1729] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1753] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1777] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1801] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1825] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1849] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1873] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1897] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1921] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1945] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1969] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1993] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [2017] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [2041] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [2065] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [2089] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [2113] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [2137] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## Levels: G0 G2
data <- data[data$Generation=="G2",]

#Remove GF and Grape
data <- data[data$Test_environment!="GF",]
data <- data[data$Test_environment!="Grape",]
data$Test_environment <- factor(data$Test_environment)
data$EggScore <- as.factor(ifelse(data$Nb_eggs<51, 1, ifelse(data$Nb_eggs<101, 2, ifelse(data$Nb_eggs<151, 3, 4))))
data$EggScoreFive <- as.factor(ifelse(data$Nb_eggs<51, 1, ifelse(data$Nb_eggs<101, 2, ifelse(data$Nb_eggs<151, 3, ifelse(data$Nb_eggs<201, 4, 5)))))
data$EggScoreSmall <- as.factor(ifelse(data$Nb_eggs<51, 1, ifelse(data$Nb_eggs<76, 2, ifelse(data$Nb_eggs<101, 3, ifelse(data$Nb_eggs<126, 4, ifelse(data$Nb_eggs<151, 5, 6))))))


## Remove WT3
data <- data[data$Population!="WT3",]
data$Original_environment <- factor(data$Original_environment)

## Remove Cherry52 which has a very high mean
#data <- data[data$Population!="Cherry52",]
data$Population <- factor(data$Population)

data$SA <- as.factor(ifelse(data$Original_environment==data$Test_environment, 1, 0))


data$Emergence_Rate <- data$Nb_adults/data$Nb_eggs
head(data)
##      Generation Experiment Original_environment Population     Date_P
## 1074         G2 Adaptation               Cherry    Cherry3 27/08/2018
## 1075         G2 Adaptation               Cherry    Cherry3 27/08/2018
## 1078         G2 Adaptation               Cherry    Cherry3 27/08/2018
## 1079         G2 Adaptation               Cherry    Cherry3 27/08/2018
## 1080         G2 Adaptation               Cherry    Cherry3 27/08/2018
## 1083         G2 Adaptation               Cherry    Cherry3 27/08/2018
##        Date_C_O Date_C_A Row Column Rack Test_environment Nb_eggs Obs_O
## 1074 28/08/2018  12/9/18  L1     C1    1           Cherry     180    LO
## 1075 28/08/2018  12/9/18  L1     C2    1       Blackberry      85    LO
## 1078 28/08/2018  12/9/18  L1     C5    1       Strawberry      67    LO
## 1079 28/08/2018  12/9/18  L2     C1    1           Cherry     106    LO
## 1080 28/08/2018  12/9/18  L2     C2    1       Blackberry     113    LO
## 1083 28/08/2018  12/9/18  L2     C5    1       Strawberry      78    LO
##      Nb_adults Obs_A EggScore EggScoreFive EggScoreSmall SA Emergence_Rate
## 1074         7    LO        4            4             6  1     0.03888889
## 1075         3    LO        2            2             3  0     0.03529412
## 1078         0    LO        2            2             2  0     0.00000000
## 1079         2    LO        3            3             4  1     0.01886792
## 1080         0    LO        3            3             4  0     0.00000000
## 1083         0    LO        2            2             3  0     0.00000000
## Compute emergence rate
data$Emergence_Rate[data$Emergence_Rate>1] <- rep(1, sum(data$Emergence_Rate>1))

2.2 Number of eggs

2.2.1 Explore data

histogram(~Nb_eggs|Original_environment*Test_environment, data=data)

tapply(data$Nb_eggs, list(data$Original_environment, data$Test_environment), mean)
##            Blackberry   Cherry Strawberry
## Blackberry   103.9720 121.8491   98.71963
## Cherry       131.7174 133.8400  103.51064
## Strawberry   119.1739 119.4894  111.56522

2.2.2 Analysis

m1 <- lm(log(Nb_eggs+1) ~ Population + Test_environment + SA + Original_environment:Test_environment, data=data[data$Generation=="G2",])
summary(m1)
## 
## Call:
## lm(formula = log(Nb_eggs + 1) ~ Population + Test_environment + 
##     SA + Original_environment:Test_environment, data = data[data$Generation == 
##     "G2", ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2545 -0.1778  0.0195  0.2192  1.0845 
## 
## Coefficients: (3 not defined because of singularities)
##                                                             Estimate Std. Error
## (Intercept)                                                4.5806543  0.1382464
## PopulationBlackberry32                                    -0.6895997  0.1067135
## PopulationBlackberry33                                     0.0419347  0.1042064
## PopulationBlackberry34                                    -0.2084134  0.1524978
## PopulationBlackberry35                                    -0.1012098  0.1033412
## PopulationBlackberry36                                    -0.0001252  0.1319869
## PopulationBlackberry37                                    -0.0843005  0.1106111
## PopulationBlackberry38                                     0.1732260  0.2076342
## PopulationBlackberry39                                     0.1099952  0.1195802
## PopulationBlackberry40                                    -0.2846811  0.1130724
## PopulationBlackberry43                                    -0.3354278  0.1578266
## PopulationBlackberry44                                     0.1972677  0.1130724
## PopulationBlackberry45                                    -0.0814644  0.2285211
## PopulationCherry103                                       -0.1405432  0.1419273
## PopulationCherry104                                       -0.8139547  0.1482000
## PopulationCherry3                                         -0.1536887  0.1750231
## PopulationCherry47                                         0.1054940  0.1228905
## PopulationCherry50                                         0.1919031  0.1635145
## PopulationCherry51                                         0.0426110  0.4444066
## PopulationCherry52                                        -0.1050116  0.2381451
## PopulationCherry6                                          0.3390668  0.1889877
## PopulationCherry7                                         -0.1152803  0.1241927
## PopulationStrawberry42                                     0.1155414  0.1196040
## PopulationStrawberry44                                     0.1533210  0.1285429
## PopulationStrawberry53                                    -0.4392990  0.1155464
## Test_environmentCherry                                     0.2938095  0.1345022
## Test_environmentStrawberry                                 0.0549272  0.1060046
## SA1                                                        0.0565484  0.1059980
## Test_environmentBlackberry:Original_environmentCherry      0.3183771  0.1384989
## Test_environmentCherry:Original_environmentCherry          0.0011227  0.1703268
## Test_environmentStrawberry:Original_environmentCherry             NA         NA
## Test_environmentBlackberry:Original_environmentStrawberry  0.2444723  0.1834154
## Test_environmentCherry:Original_environmentStrawberry             NA         NA
## Test_environmentStrawberry:Original_environmentStrawberry         NA         NA
##                                                           t value Pr(>|t|)    
## (Intercept)                                                33.134  < 2e-16 ***
## PopulationBlackberry32                                     -6.462 2.21e-10 ***
## PopulationBlackberry33                                      0.402 0.687526    
## PopulationBlackberry34                                     -1.367 0.172268    
## PopulationBlackberry35                                     -0.979 0.327809    
## PopulationBlackberry36                                     -0.001 0.999244    
## PopulationBlackberry37                                     -0.762 0.446294    
## PopulationBlackberry38                                      0.834 0.404470    
## PopulationBlackberry39                                      0.920 0.358043    
## PopulationBlackberry40                                     -2.518 0.012085 *  
## PopulationBlackberry43                                     -2.125 0.033991 *  
## PopulationBlackberry44                                      1.745 0.081590 .  
## PopulationBlackberry45                                     -0.356 0.721609    
## PopulationCherry103                                        -0.990 0.322473    
## PopulationCherry104                                        -5.492 5.98e-08 ***
## PopulationCherry3                                          -0.878 0.380256    
## PopulationCherry47                                          0.858 0.391010    
## PopulationCherry50                                          1.174 0.241038    
## PopulationCherry51                                          0.096 0.923647    
## PopulationCherry52                                         -0.441 0.659411    
## PopulationCherry6                                           1.794 0.073323 .  
## PopulationCherry7                                          -0.928 0.353677    
## PopulationStrawberry42                                      0.966 0.334436    
## PopulationStrawberry44                                      1.193 0.233458    
## PopulationStrawberry53                                     -3.802 0.000159 ***
## Test_environmentCherry                                      2.184 0.029337 *  
## Test_environmentStrawberry                                  0.518 0.604548    
## SA1                                                         0.533 0.593905    
## Test_environmentBlackberry:Original_environmentCherry       2.299 0.021878 *  
## Test_environmentCherry:Original_environmentCherry           0.007 0.994743    
## Test_environmentStrawberry:Original_environmentCherry          NA       NA    
## Test_environmentBlackberry:Original_environmentStrawberry   1.333 0.183100    
## Test_environmentCherry:Original_environmentStrawberry          NA       NA    
## Test_environmentStrawberry:Original_environmentStrawberry      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4262 on 571 degrees of freedom
## Multiple R-squared:  0.3375, Adjusted R-squared:  0.3027 
## F-statistic: 9.697 on 30 and 571 DF,  p-value: < 2.2e-16
anova(m1)
## Analysis of Variance Table
## 
## Response: log(Nb_eggs + 1)
##                                        Df  Sum Sq Mean Sq F value    Pr(>F)    
## Population                             24  45.525 1.89686 10.4403 < 2.2e-16 ***
## Test_environment                        2   5.834 2.91691 16.0546 1.647e-07 ***
## SA                                      1   0.404 0.40421  2.2248    0.1364    
## Test_environment:Original_environment   3   1.093 0.36447  2.0060    0.1120    
## Residuals                             571 103.743 0.18169                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## F test for SA
(Fratio_int = (anova(m1)[3,2]/anova(m1)[4,2])/(1/anova(m1)[4,1]))
## [1] 1.10903
(pvalue_int = 1 - pf(Fratio_int, 1, anova(m1)[4,1])) #the correct test (see equation D7 in Appendix D of the paper)
## [1] 0.3696212
## Compute R2 = MS Interaction model without SA - MS Interaction model with SA / MS Interaction model without SA
(rsq2 <- 1-anova(m1)[4, 3]/((anova(m1)[3, 2]+anova(m1)[4, 2])/(anova(m1)[3, 1]+anova(m1)[4, 1])))
## [1] 0.02653419

2.3 Number of adults

2.3.1 Explore data with blue

m0 <- lm(Nb_adults~Original_environment*Test_environment, data=data[data$Generation=="G2",])

summary(m0)
## 
## Call:
## lm(formula = Nb_adults ~ Original_environment * Test_environment, 
##     data = data[data$Generation == "G2", ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -27.318 -11.287  -2.906   8.834  85.087 
## 
## Coefficients:
##                                                           Estimate Std. Error
## (Intercept)                                                27.3178     1.5742
## Original_environmentCherry                                 -4.4047     2.8710
## Original_environmentStrawberry                             -1.6656     2.8710
## Test_environmentCherry                                     -4.1291     2.2315
## Test_environmentStrawberry                                -14.1589     2.2263
## Original_environmentCherry:Test_environmentCherry           4.5360     4.0059
## Original_environmentStrawberry:Test_environmentCherry       0.8386     4.0479
## Original_environmentCherry:Test_environmentStrawberry       2.7777     4.0450
## Original_environmentStrawberry:Test_environmentStrawberry   4.5502     4.0602
##                                                           t value Pr(>|t|)    
## (Intercept)                                                17.353  < 2e-16 ***
## Original_environmentCherry                                 -1.534   0.1255    
## Original_environmentStrawberry                             -0.580   0.5620    
## Test_environmentCherry                                     -1.850   0.0648 .  
## Test_environmentStrawberry                                 -6.360 4.03e-10 ***
## Original_environmentCherry:Test_environmentCherry           1.132   0.2579    
## Original_environmentStrawberry:Test_environmentCherry       0.207   0.8359    
## Original_environmentCherry:Test_environmentStrawberry       0.687   0.4925    
## Original_environmentStrawberry:Test_environmentStrawberry   1.121   0.2629    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 16.28 on 593 degrees of freedom
## Multiple R-squared:  0.1037, Adjusted R-squared:  0.09163 
## F-statistic: 8.578 on 8 and 593 DF,  p-value: 4.287e-11
tapply(data$Nb_adults[data$Generation=="G2"], list(data$Original_environment[data$Generation=="G2"], data$Test_environment[data$Generation=="G2"]), mean)
##            Blackberry   Cherry Strawberry
## Blackberry   27.31776 23.18868   13.15888
## Cherry       22.91304 23.32000   11.53191
## Strawberry   25.65217 22.36170   16.04348
tapply(data$Nb_adults[data$Generation=="G2"], list(data$Original_environment[data$Generation=="G2"], data$Test_environment[data$Generation=="G2"]), var)
##            Blackberry   Cherry Strawberry
## Blackberry   402.3698 205.4688   143.9085
## Cherry       506.6589 412.7935   174.1240
## Strawberry   337.0763 133.4968   120.3092
range(data$Nb_adults)
## [1]   0 108
## Check number of eggs and adults
tapply(data$Nb_eggs[data$Generation=="G2"], list(data$Population[data$Generation=="G2"], data$Test_environment[data$Generation=="G2"]), mean)
##              Blackberry    Cherry Strawberry
## Blackberry31  106.44444 130.33333  118.88889
## Blackberry32   61.23077  84.46154   74.00000
## Blackberry33  125.13333 141.28571   95.66667
## Blackberry34  105.75000  95.66667   83.50000
## Blackberry35   93.60000 123.73333   96.31250
## Blackberry36  118.00000 144.00000   92.83333
## Blackberry37  108.09091 114.45455  101.27273
## Blackberry38  125.50000 141.00000  143.00000
## Blackberry39  133.75000 146.00000  117.12500
## Blackberry40   88.30000  92.10000   80.40000
## Blackberry43   69.33333  90.00000   86.33333
## Blackberry44  129.20000 155.70000  129.70000
## Blackberry45   92.00000 113.00000  110.00000
## Cherry103     129.66667 118.00000   89.83333
## Cherry104      59.80000  61.83333   55.20000
## Cherry3       109.66667 143.66667   78.00000
## Cherry47      142.00000 162.71429  128.23077
## Cherry50      148.33333 185.00000  126.00000
## Cherry51      139.00000        NA         NA
## Cherry52      129.50000  96.00000  103.00000
## Cherry6       226.50000 203.50000  124.00000
## Cherry7       137.75000 118.69231   97.50000
## Strawberry42  140.46667 131.81250  146.66667
## Strawberry44  150.60000 142.30000  139.70000
## Strawberry53   89.00000  99.23810   73.09524
tapply(data$Nb_adults[data$Generation=="G2"], list(data$Population[data$Generation=="G2"], data$Test_environment[data$Generation=="G2"]), mean)
##              Blackberry    Cherry Strawberry
## Blackberry31  38.444444 36.777778 19.2222222
## Blackberry32   1.923077  6.307692  0.6923077
## Blackberry33  30.800000 34.000000 16.6666667
## Blackberry34  13.500000 17.666667  8.0000000
## Blackberry35  14.533333 23.533333  5.2500000
## Blackberry36  28.800000 23.166667 14.0000000
## Blackberry37  34.727273 24.090909  8.2727273
## Blackberry38  15.500000  7.000000  1.0000000
## Blackberry39  28.375000 26.750000 19.0000000
## Blackberry40  38.700000 17.300000 23.9000000
## Blackberry43  45.000000 13.750000 11.6666667
## Blackberry44  48.900000 27.300000 25.6000000
## Blackberry45  11.500000 30.000000  2.0000000
## Cherry103     35.666667 41.857143 20.8333333
## Cherry104      8.200000 10.833333  3.0000000
## Cherry3        1.000000  3.000000  0.0000000
## Cherry47      15.833333 18.571429  5.9230769
## Cherry50      63.000000 67.750000 46.0000000
## Cherry51      78.000000        NA         NA
## Cherry52      65.500000 39.000000 18.0000000
## Cherry6        5.500000 10.000000  1.3333333
## Cherry7       16.416667 16.076923  9.9166667
## Strawberry42  45.733333 35.000000 26.3333333
## Strawberry44  13.100000 13.400000  9.9000000
## Strawberry53  17.285714 17.000000 11.6190476
## Check for the presence of negative correlations
m1 <- lm(log(Nb_adults+1) ~ Population + Test_environment, data=data[data$Generation=="G2",])
data$resid <- residuals(m1)

meanbypopbytestenv <- as.data.frame(tapply(data$resid[data$Generation=="G2"], list(data$Population[data$Generation=="G2"], data$Test_environment[data$Generation=="G2"]), mean))


## Cherry ~ Blackberry
plot(meanbypopbytestenv$Blackberry, meanbypopbytestenv$Cherry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), xlab="Blackberry", ylab="Cherry", pch=16)
legend("topright", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)

## Strawberry ~ Cherry
plot(meanbypopbytestenv$Cherry, meanbypopbytestenv$Strawberry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), xlab="Cherry", ylab="Strawberry", pch=16)
legend("topright", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)

## Strawberry ~ Blackberry
plot(meanbypopbytestenv$Blackberry~ meanbypopbytestenv$Strawberry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), xlab="Blackberry", ylab="Strawberry", pch=16)
legend("topright", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)

2.3.2 Explore data with BLUPS

m0 <- lmer(log(Nb_adults+1)~ Population + Test_environment + (0 + Test_environment | Population), data = data[data$Generation=="G2",])
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues
## Correlation between fixed effects
plot(coef(m1), fixef(m0))

summary(m0)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## log(Nb_adults + 1) ~ Population + Test_environment + (0 + Test_environment |  
##     Population)
##    Data: data[data$Generation == "G2", ]
## 
## REML criterion at convergence: 1264.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8381 -0.4992  0.1024  0.5908  2.6516 
## 
## Random effects:
##  Groups     Name                       Variance  Std.Dev. Corr       
##  Population Test_environmentBlackberry 9.254e-07 0.000962            
##             Test_environmentCherry     8.208e-02 0.286500  0.59      
##             Test_environmentStrawberry 4.287e-02 0.207058  0.44 -0.13
##  Residual                              4.404e-01 0.663613            
## Number of obs: 602, groups:  Population, 25
## 
## Fixed effects:
##                            Estimate Std. Error t value
## (Intercept)                 3.55606    0.16116  22.066
## PopulationBlackberry32     -2.50179    0.20907 -11.966
## PopulationBlackberry33     -0.17462    0.20469  -0.853
## PopulationBlackberry34     -1.17436    0.27159  -4.324
## PopulationBlackberry35     -1.04239    0.20424  -5.104
## PopulationBlackberry36     -0.32019    0.24803  -1.291
## PopulationBlackberry37     -0.52201    0.21491  -2.429
## PopulationBlackberry38     -1.36372    0.35317  -3.861
## PopulationBlackberry39     -0.34622    0.22792  -1.519
## PopulationBlackberry40     -0.09950    0.21853  -0.455
## PopulationBlackberry43     -0.48439    0.28519  -1.698
## PopulationBlackberry44      0.16978    0.21853   0.777
## PopulationBlackberry45     -1.09654    0.37593  -2.917
## PopulationCherry103         0.08807    0.24101   0.365
## PopulationCherry104        -1.51454    0.25101  -6.034
## PopulationCherry3          -2.86774    0.29040  -9.875
## PopulationCherry47         -1.07034    0.21085  -5.076
## PopulationCherry50          0.68579    0.27800   2.467
## PopulationCherry51          0.81338    0.68290   1.191
## PopulationCherry52          0.37570    0.37593   0.999
## PopulationCherry6          -1.98132    0.31557  -6.279
## PopulationCherry7          -0.82275    0.21155  -3.889
## PopulationStrawberry42      0.23210    0.20440   1.136
## PopulationStrawberry44     -0.93199    0.21853  -4.265
## PopulationStrawberry53     -0.75587    0.19532  -3.870
## Test_environmentCherry      0.04809    0.09288   0.518
## Test_environmentStrawberry -0.70086    0.08215  -8.531
## 
## Correlation matrix not shown by default, as p = 27 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
## optimizer (nloptwrap) convergence code: 0 (OK)
## unable to evaluate scaled gradient
## Model failed to converge: degenerate  Hessian with 1 negative eigenvalues
BLUP <- ranef(m0)$Population

pairs(BLUP)

## Negative correlation between blups and blues, very weird, no?
plot(BLUP$Test_environmentBlackberry, meanbypopbytestenv$Blackberry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), pch=16)
legend("topright", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)

## Positive correlation betwen blups and blues, as expected
plot(BLUP$Test_environmentCherry, meanbypopbytestenv$Cherry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), pch=16)
legend("topleft", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)

## Positive correlation betwen blups and blues, as expected
plot(BLUP$Test_environmentStrawberry, meanbypopbytestenv$Strawberry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), pch=16)
legend("topright", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)

2.3.3 Analysis

## Null model
m0 <- lm(log(Nb_adults+1) ~ Population + Test_environment + Original_environment:Test_environment, data=data[data$Generation=="G2",])

## SA model
m1 <- lm(log(Nb_adults+1) ~ Population + Test_environment + SA + Original_environment:Test_environment, data=data[data$Generation=="G2",])
summary(m1)
## 
## Call:
## lm(formula = log(Nb_adults + 1) ~ Population + Test_environment + 
##     SA + Original_environment:Test_environment, data = data[data$Generation == 
##     "G2", ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.50160 -0.35356  0.05873  0.42543  1.98747 
## 
## Coefficients: (3 not defined because of singularities)
##                                                           Estimate Std. Error
## (Intercept)                                                3.20150    0.22142
## PopulationBlackberry32                                    -2.35665    0.17091
## PopulationBlackberry33                                    -0.12012    0.16690
## PopulationBlackberry34                                    -1.08199    0.24424
## PopulationBlackberry35                                    -0.98458    0.16551
## PopulationBlackberry36                                    -0.32673    0.21139
## PopulationBlackberry37                                    -0.61653    0.17716
## PopulationBlackberry38                                    -1.39482    0.33255
## PopulationBlackberry39                                    -0.31957    0.19152
## PopulationBlackberry40                                    -0.16801    0.18110
## PopulationBlackberry43                                    -0.57342    0.25278
## PopulationBlackberry44                                     0.11078    0.18110
## PopulationBlackberry45                                    -1.06627    0.36600
## PopulationCherry103                                        0.11034    0.22731
## PopulationCherry104                                       -1.47775    0.23736
## PopulationCherry3                                         -2.81107    0.28032
## PopulationCherry47                                        -1.04591    0.19682
## PopulationCherry50                                         0.71795    0.26189
## PopulationCherry51                                         0.88670    0.71177
## PopulationCherry52                                         0.39346    0.38142
## PopulationCherry6                                         -1.92021    0.30268
## PopulationCherry7                                         -0.78567    0.19891
## PopulationStrawberry42                                     0.08779    0.19156
## PopulationStrawberry44                                    -1.01054    0.20588
## PopulationStrawberry53                                    -0.82355    0.18506
## Test_environmentCherry                                     0.38655    0.21542
## Test_environmentStrawberry                                -0.42287    0.16978
## SA1                                                        0.38209    0.16977
## Test_environmentBlackberry:Original_environmentCherry      0.28124    0.22182
## Test_environmentCherry:Original_environmentCherry         -0.27686    0.27280
## Test_environmentStrawberry:Original_environmentCherry           NA         NA
## Test_environmentBlackberry:Original_environmentStrawberry  0.39178    0.29376
## Test_environmentCherry:Original_environmentStrawberry           NA         NA
## Test_environmentStrawberry:Original_environmentStrawberry       NA         NA
##                                                           t value Pr(>|t|)    
## (Intercept)                                                14.459  < 2e-16 ***
## PopulationBlackberry32                                    -13.789  < 2e-16 ***
## PopulationBlackberry33                                     -0.720 0.471991    
## PopulationBlackberry34                                     -4.430 1.13e-05 ***
## PopulationBlackberry35                                     -5.949 4.72e-09 ***
## PopulationBlackberry36                                     -1.546 0.122748    
## PopulationBlackberry37                                     -3.480 0.000539 ***
## PopulationBlackberry38                                     -4.194 3.17e-05 ***
## PopulationBlackberry39                                     -1.669 0.095747 .  
## PopulationBlackberry40                                     -0.928 0.353952    
## PopulationBlackberry43                                     -2.268 0.023673 *  
## PopulationBlackberry44                                      0.612 0.540970    
## PopulationBlackberry45                                     -2.913 0.003716 ** 
## PopulationCherry103                                         0.485 0.627572    
## PopulationCherry104                                        -6.226 9.29e-10 ***
## PopulationCherry3                                         -10.028  < 2e-16 ***
## PopulationCherry47                                         -5.314 1.54e-07 ***
## PopulationCherry50                                          2.741 0.006308 ** 
## PopulationCherry51                                          1.246 0.213359    
## PopulationCherry52                                          1.032 0.302711    
## PopulationCherry6                                          -6.344 4.56e-10 ***
## PopulationCherry7                                          -3.950 8.80e-05 ***
## PopulationStrawberry42                                      0.458 0.646932    
## PopulationStrawberry44                                     -4.909 1.20e-06 ***
## PopulationStrawberry53                                     -4.450 1.03e-05 ***
## Test_environmentCherry                                      1.794 0.073279 .  
## Test_environmentStrawberry                                 -2.491 0.013032 *  
## SA1                                                         2.251 0.024786 *  
## Test_environmentBlackberry:Original_environmentCherry       1.268 0.205355    
## Test_environmentCherry:Original_environmentCherry          -1.015 0.310591    
## Test_environmentStrawberry:Original_environmentCherry          NA       NA    
## Test_environmentBlackberry:Original_environmentStrawberry   1.334 0.182843    
## Test_environmentCherry:Original_environmentStrawberry          NA       NA    
## Test_environmentStrawberry:Original_environmentStrawberry      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6827 on 571 degrees of freedom
## Multiple R-squared:  0.5989, Adjusted R-squared:  0.5778 
## F-statistic: 28.42 on 30 and 571 DF,  p-value: < 2.2e-16
anova(m1)
## Analysis of Variance Table
## 
## Response: log(Nb_adults + 1)
##                                        Df Sum Sq Mean Sq F value    Pr(>F)    
## Population                             24 323.52  13.480 28.9240 < 2.2e-16 ***
## Test_environment                        2  69.66  34.828 74.7289 < 2.2e-16 ***
## SA                                      1   3.21   3.215  6.8980  0.008861 ** 
## Test_environment:Original_environment   3   0.99   0.328  0.7046  0.549525    
## Residuals                             571 266.12   0.466                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## F test for SA
(Fratio_int = (anova(m1)[3,2]/anova(m1)[4,2])/(1/anova(m1)[4,1]))
## [1] 9.789562
(pvalue_int = 1 - pf(Fratio_int, 1, anova(m1)[4,1])) #the correct test (see equation D7 in Appendix D of the paper)
## [1] 0.05211294
## Compute R2 = MS Interaction model without SA - MS Interaction model with SA / MS Interaction model without SA
## rsquare with two models
rsq2 <- (anova(m0)[3, 3]-anova(m1)[4, 3])/anova(m0)[3, 3]
## rsquare with single model
(rsq2 <- 1-anova(m1)[4, 3]/((anova(m1)[3, 2]+anova(m1)[4, 2])/(anova(m1)[3, 1]+anova(m1)[4, 1])))
## [1] 0.687245

2.3.4 Power Analysis

m2 <- lm(log(Nb_adults+1) ~ -1 + Population + Test_environment + Original_environment:Test_environment, data=data[data$Generation=="G2",])
coef(m2)

simdata <- function(seed=1, nrep = 2, npop = 3, SAincrease=0){
  set.seed(seed)
  ## Sample populations
  PopBlackberry <- sort(sample(levels(data$Population)[grep("Blackberry", levels(data$Population))], size=npop, replace = TRUE))
  PopCherry <- sort(sample(levels(data$Population)[grep("Cherry", levels(data$Population))], size=npop, replace = TRUE))
  PopStrawberry <- sort(sample(levels(data$Population)[grep("Strawberry", levels(data$Population))], size=npop, replace = TRUE))
  
  ## Create dataset
  newdata <- expand.grid(Pop_new_name=1:(3*npop), Test_environment=levels(data$Test_environment), Rep=1:nrep)
  newdata$Original_environment <- rep(c("Blackberry", "Cherry", "Strawberry"), each=npop)
  
  ## Add original names to get the right coef in the model
  newdata$Pop_new_name <- as.factor(newdata$Pop_new_name)
  newdata$Population <- newdata$Pop_new_name
  levels(newdata$Population) <- c(PopBlackberry, PopCherry, PopStrawberry)
  
  ## Get design matrix
  mat <- model.matrix(~-1 + Population + Test_environment + Original_environment:Test_environment, data=newdata)
  
  ## Check that the order of the coefficients is the same
  data.frame(colnames(mat), names(coef(m2))[names(coef(m2))%in%colnames(mat)])
  
  coefsim <- coef(m2)[names(coef(m2))%in%colnames(mat)]
  coefsim <- ifelse(is.na(coefsim), 0, coefsim)
  
  ## Increase fitness in sympatric environment
  coefsim[grepl("PopulationBlackberry|Test_environmentCherry:Original_environmentCherry|Test_environmentStrawberry:Original_environmentStrawberry", names(coefsim))] <- coefsim[grepl("PopulationBlackberry|Test_environmentCherry:Original_environmentCherry|Test_environmentStrawberry:Original_environmentStrawberry", names(coefsim))]+SAincrease
  
  newdata$resp <- mat%*%coefsim + rnorm(n=nrow(newdata), mean = 0, sd = sigma(m2))
  
  m3 <- lm(resp ~ -1+Population + Test_environment + Original_environment:Test_environment, data=newdata)
  data.frame(coefsim, coef(m3))
  
  newdata$SA <- as.factor(ifelse(newdata$Original_environment==newdata$Test_environment, 1, 0))
  
  m4 <- lm(resp ~ Pop_new_name + Test_environment + SA + Original_environment:Test_environment, data=newdata)
  summary(m4)
  anova(m4)
    
  ## F test for SA
  (Fratio_int = (anova(m4)[3,2]/anova(m4)[4,2])/(1/anova(m4)[4,1]))
  (pvalue_int = 1 - pf(Fratio_int, 1, anova(m4)[4,1])) #the correct test (see equation D7 in Appendix D of the paper)
  ## rsquare with single model
  (rsq2 <- 1-anova(m4)[4, 3]/((anova(m4)[3, 2]+anova(m4)[4, 2])/(anova(m4)[3, 1]+anova(m4)[4, 1])))
  #print(rsq2)
  return(pvalue_int)
}

simdata(seed=2, nrep = 30, npop = 3)
simdata(seed=2, nrep = 30, npop = 3)

simpval <- sapply(1:500, simdata, nrep = 10, npop = 3)
mean(simpval)
## Compute power
mean(simpval<0.05)
hist(simpval)

simpval10rep <- sapply(1:500, simdata, nrep = 10, npop = 3)
simpval15rep <- sapply(1:500, simdata, nrep = 15, npop = 3)
simpval20rep <- sapply(1:500, simdata, nrep = 20, npop = 3)
simpval25rep <- sapply(1:500, simdata, nrep = 25, npop = 3)
simpval30rep <- sapply(1:500, simdata, nrep = 30, npop = 3)
simpval35rep <- sapply(1:500, simdata, nrep = 35, npop = 3)
simpval40rep <- sapply(1:500, simdata, nrep = 40, npop = 3)
simpval50rep <- sapply(1:500, simdata, nrep = 50, npop = 3)
simpval60rep <- sapply(1:500, simdata, nrep = 60, npop = 3)
simpval80rep <- sapply(1:500, simdata, nrep = 80, npop = 3)
simpval100rep <- sapply(1:500, simdata, nrep = 100, npop = 3)

## Compute power
val <- data.frame(simpval10rep, simpval15rep, simpval20rep, simpval25rep, simpval30rep, simpval35rep, simpval40rep, simpval50rep, simpval60rep, simpval80rep, simpval100rep)

computepower <- function(x){
  return(mean(x<0.05))
}

statpower <- apply(val, 2, computepower)

plot(c(10, 15, 20, 25, 30, 35, 40, 50, 60, 80, 100), statpower, las=1, xlab="Number of replicates", ylab="Statistical Power")

2.4 Number of adults~ Number of eggs

2.4.1 Explore data

AvNb_adults <- tapply(data$Nb_adults, list(data$EggScoreFive, data$Original_environment, data$Test_environment), mean)

AvNb_adults[, , "Blackberry"][, "Cherry"]/AvNb_adults[, , "Blackberry"][, "Blackberry"]
##         1         2         3         4         5 
## 0.2796610 0.7280423 0.7701180 1.0788644 0.1762673
AvNb_adults[, , "Cherry"][, "Blackberry"]/AvNb_adults[, , "Cherry"][, "Cherry"]
##         1         2         3         4         5 
## 0.0000000 0.8869732 1.1581395 0.8547641 1.5423729
xyplot(Nb_adults~Nb_eggs|Original_environment*Test_environment, data=data)

m0 <- lm(log(Nb_adults+1)~log(Nb_eggs+1)+Test_environment+Original_environment, data=data)
m1 <- lm(log(Nb_adults+1)~log(Nb_eggs+1)*Test_environment+Original_environment, data=data)
m2 <- lm(log(Nb_adults+1)~log(Nb_eggs+1)*Original_environment+Test_environment, data=data)

m3 <- lm(log(Nb_adults+1)~Nb_eggs+Test_environment+Original_environment, data=data)

summary(m1)
## 
## Call:
## lm(formula = log(Nb_adults + 1) ~ log(Nb_eggs + 1) * Test_environment + 
##     Original_environment, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3193 -0.4821  0.1630  0.6116  1.7523 
## 
## Coefficients:
##                                             Estimate Std. Error t value
## (Intercept)                                 -0.69352    0.50780  -1.366
## log(Nb_eggs + 1)                             0.78227    0.10945   7.147
## Test_environmentCherry                      -0.69176    1.01112  -0.684
## Test_environmentStrawberry                  -0.57547    0.74678  -0.771
## Original_environmentCherry                  -0.29381    0.09165  -3.206
## Original_environmentStrawberry               0.15817    0.09208   1.718
## log(Nb_eggs + 1):Test_environmentCherry      0.13184    0.21267   0.620
## log(Nb_eggs + 1):Test_environmentStrawberry -0.01274    0.16191  -0.079
##                                             Pr(>|t|)    
## (Intercept)                                  0.17254    
## log(Nb_eggs + 1)                            2.61e-12 ***
## Test_environmentCherry                       0.49415    
## Test_environmentStrawberry                   0.44125    
## Original_environmentCherry                   0.00142 ** 
## Original_environmentStrawberry               0.08635 .  
## log(Nb_eggs + 1):Test_environmentCherry      0.53556    
## log(Nb_eggs + 1):Test_environmentStrawberry  0.93733    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.904 on 594 degrees of freedom
## Multiple R-squared:  0.2683, Adjusted R-squared:  0.2597 
## F-statistic: 31.12 on 7 and 594 DF,  p-value: < 2.2e-16
##Cl: density dependent larval survival rate does not seem to depend on the environment
anova(m1)
## Analysis of Variance Table
## 
## Response: log(Nb_adults + 1)
##                                    Df Sum Sq Mean Sq  F value    Pr(>F)    
## log(Nb_eggs + 1)                    1 114.92 114.919 140.6112 < 2.2e-16 ***
## Test_environment                    2  47.80  23.898  29.2403 7.727e-13 ***
## Original_environment                2  14.92   7.462   9.1301 0.0001243 ***
## log(Nb_eggs + 1):Test_environment   2   0.39   0.197   0.2406 0.7862571    
## Residuals                         594 485.47   0.817                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(m2)
## Analysis of Variance Table
## 
## Response: log(Nb_adults + 1)
##                                        Df Sum Sq Mean Sq  F value    Pr(>F)    
## log(Nb_eggs + 1)                        1 114.92 114.919 141.5715 < 2.2e-16 ***
## Original_environment                    2  15.32   7.661   9.4372 9.231e-05 ***
## Test_environment                        2  47.40  23.699  29.1953 8.050e-13 ***
## log(Nb_eggs + 1):Original_environment   2   3.69   1.843   2.2706    0.1041    
## Residuals                             594 482.17   0.812                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(MuMIn)
model.sel(m0, m1, m2, m3)
## Model selection table 
##      (Int) log(Nb_egg+1) Org_env Tst_env log(Nb_egg+1):Tst_env
## m2 -1.3480        0.9262       +       +                      
## m0 -0.7699        0.7988       +       +                      
## m1 -0.6935        0.7823       +       +                     +
## m3  1.9780                     +       +                      
##    log(Nb_egg+1):Org_env   Nb_egg df   logLik   AICc delta weight
## m2                     +           9 -787.393 1593.1  0.00  0.521
## m0                                 7 -789.685 1593.6  0.47  0.412
## m1                                 9 -789.441 1597.2  4.10  0.067
## m3                       0.008389  7 -801.139 1616.5 23.38  0.000
## Models ranked by AICc(x)

2.4.2 Analysis

m1 <- lm(log(Nb_adults+1) ~ Population + Test_environment + SA + Original_environment:Test_environment+log(Nb_eggs+1), data=data[data$Generation=="G2",])
summary(m1)
## 
## Call:
## lm(formula = log(Nb_adults + 1) ~ Population + Test_environment + 
##     SA + Original_environment:Test_environment + log(Nb_eggs + 
##     1), data = data[data$Generation == "G2", ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.27701 -0.34196  0.05666  0.38365  1.72524 
## 
## Coefficients: (3 not defined because of singularities)
##                                                           Estimate Std. Error
## (Intercept)                                                0.86228    0.35909
## PopulationBlackberry32                                    -2.00449    0.16796
## PopulationBlackberry33                                    -0.14154    0.15835
## PopulationBlackberry34                                    -0.97556    0.23208
## PopulationBlackberry35                                    -0.93289    0.15714
## PopulationBlackberry36                                    -0.32667    0.20053
## PopulationBlackberry37                                    -0.57348    0.16814
## PopulationBlackberry38                                    -1.48328    0.31566
## PopulationBlackberry39                                    -0.37574    0.18182
## PopulationBlackberry40                                    -0.02263    0.17275
## PopulationBlackberry43                                    -0.40212    0.24074
## PopulationBlackberry44                                     0.01004    0.17225
## PopulationBlackberry45                                    -1.02467    0.34724
## PopulationCherry103                                        0.18211    0.21582
## PopulationCherry104                                       -1.06208    0.23104
## PopulationCherry3                                         -2.73259    0.26610
## PopulationCherry47                                        -1.09979    0.18683
## PopulationCherry50                                         0.61995    0.24873
## PopulationCherry51                                         0.86494    0.67521
## PopulationCherry52                                         0.44708    0.36189
## PopulationCherry6                                         -2.09337    0.28795
## PopulationCherry7                                         -0.72679    0.18883
## PopulationStrawberry42                                     0.02878    0.18187
## PopulationStrawberry44                                    -1.08884    0.19554
## PopulationStrawberry53                                    -0.59921    0.17776
## Test_environmentCherry                                     0.23651    0.20521
## Test_environmentStrawberry                                -0.45092    0.16110
## SA1                                                        0.35321    0.16109
## log(Nb_eggs + 1)                                           0.51067    0.06358
## Test_environmentBlackberry:Original_environmentCherry      0.11866    0.21140
## Test_environmentCherry:Original_environmentCherry         -0.27743    0.25879
## Test_environmentStrawberry:Original_environmentCherry           NA         NA
## Test_environmentBlackberry:Original_environmentStrawberry  0.26693    0.27910
## Test_environmentCherry:Original_environmentStrawberry           NA         NA
## Test_environmentStrawberry:Original_environmentStrawberry       NA         NA
##                                                           t value Pr(>|t|)    
## (Intercept)                                                 2.401 0.016656 *  
## PopulationBlackberry32                                    -11.934  < 2e-16 ***
## PopulationBlackberry33                                     -0.894 0.371793    
## PopulationBlackberry34                                     -4.204 3.05e-05 ***
## PopulationBlackberry35                                     -5.937 5.06e-09 ***
## PopulationBlackberry36                                     -1.629 0.103865    
## PopulationBlackberry37                                     -3.411 0.000694 ***
## PopulationBlackberry38                                     -4.699 3.28e-06 ***
## PopulationBlackberry39                                     -2.067 0.039225 *  
## PopulationBlackberry40                                     -0.131 0.895838    
## PopulationBlackberry43                                     -1.670 0.095398 .  
## PopulationBlackberry44                                      0.058 0.953538    
## PopulationBlackberry45                                     -2.951 0.003299 ** 
## PopulationCherry103                                         0.844 0.399132    
## PopulationCherry104                                        -4.597 5.28e-06 ***
## PopulationCherry3                                         -10.269  < 2e-16 ***
## PopulationCherry47                                         -5.886 6.74e-09 ***
## PopulationCherry50                                          2.492 0.012970 *  
## PopulationCherry51                                          1.281 0.200719    
## PopulationCherry52                                          1.235 0.217181    
## PopulationCherry6                                          -7.270 1.19e-12 ***
## PopulationCherry7                                          -3.849 0.000132 ***
## PopulationStrawberry42                                      0.158 0.874310    
## PopulationStrawberry44                                     -5.568 3.97e-08 ***
## PopulationStrawberry53                                     -3.371 0.000800 ***
## Test_environmentCherry                                      1.153 0.249586    
## Test_environmentStrawberry                                 -2.799 0.005299 ** 
## SA1                                                         2.193 0.028734 *  
## log(Nb_eggs + 1)                                            8.032 5.54e-15 ***
## Test_environmentBlackberry:Original_environmentCherry       0.561 0.574817    
## Test_environmentCherry:Original_environmentCherry          -1.072 0.284151    
## Test_environmentStrawberry:Original_environmentCherry          NA       NA    
## Test_environmentBlackberry:Original_environmentStrawberry   0.956 0.339279    
## Test_environmentCherry:Original_environmentStrawberry          NA       NA    
## Test_environmentStrawberry:Original_environmentStrawberry      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6476 on 570 degrees of freedom
## Multiple R-squared:  0.6397, Adjusted R-squared:  0.6201 
## F-statistic: 32.64 on 31 and 570 DF,  p-value: < 2.2e-16
anova(m1)
## Analysis of Variance Table
## 
## Response: log(Nb_adults + 1)
##                                        Df Sum Sq Mean Sq F value    Pr(>F)    
## Population                             24 323.52  13.480 32.1410 < 2.2e-16 ***
## Test_environment                        2  69.66  34.828 83.0404 < 2.2e-16 ***
## SA                                      1   3.21   3.215  7.6652  0.005813 ** 
## log(Nb_eggs + 1)                        1  27.52  27.522 65.6206 3.338e-15 ***
## Test_environment:Original_environment   3   0.52   0.173  0.4120  0.744439    
## Residuals                             570 239.06   0.419                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## F test for SA
(Fratio_int = (anova(m1)[3,2]/anova(m1)[5,2])/(1/anova(m1)[5,1]))
## [1] 18.60447
(pvalue_int = 1 - pf(Fratio_int, 1, anova(m1)[5,1])) #the correct test (see equation D7 in Appendix D of the paper)
## [1] 0.02295088
## Compute R2 = MS Interaction model without SA - MS Interaction model with SA / MS Interaction model without SA
(rsq2 <- 1-anova(m1)[5, 3]/((anova(m1)[3, 2]+anova(m1)[5, 2])/(anova(m1)[3, 1]+anova(m1)[5, 1])))
## [1] 0.8148531
m2 <- lm(log(Nb_adults+1) ~ Population + Test_environment + SA + Original_environment:Test_environment+EggScore, data=data[data$Generation=="G2",])
summary(m2)
## 
## Call:
## lm(formula = log(Nb_adults + 1) ~ Population + Test_environment + 
##     SA + Original_environment:Test_environment + EggScore, data = data[data$Generation == 
##     "G2", ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3975 -0.3259  0.0477  0.3925  1.6645 
## 
## Coefficients: (3 not defined because of singularities)
##                                                           Estimate Std. Error
## (Intercept)                                                2.42964    0.24025
## PopulationBlackberry32                                    -2.07722    0.16498
## PopulationBlackberry33                                    -0.17843    0.15808
## PopulationBlackberry34                                    -1.08307    0.23151
## PopulationBlackberry35                                    -0.96881    0.15722
## PopulationBlackberry36                                    -0.37498    0.20035
## PopulationBlackberry37                                    -0.59598    0.16811
## PopulationBlackberry38                                    -1.52552    0.31662
## PopulationBlackberry39                                    -0.40816    0.18147
## PopulationBlackberry40                                    -0.08035    0.17306
## PopulationBlackberry43                                    -0.39924    0.24076
## PopulationBlackberry44                                    -0.02872    0.17257
## PopulationBlackberry45                                    -1.14581    0.34702
## PopulationCherry103                                        0.12017    0.21574
## PopulationCherry104                                       -1.10701    0.22934
## PopulationCherry3                                         -2.79795    0.26564
## PopulationCherry47                                        -1.14553    0.18664
## PopulationCherry50                                         0.56136    0.24893
## PopulationCherry51                                         0.81588    0.67410
## PopulationCherry52                                         0.37138    0.36179
## PopulationCherry6                                         -2.06312    0.28701
## PopulationCherry7                                         -0.76473    0.18843
## PopulationStrawberry42                                    -0.01311    0.18156
## PopulationStrawberry44                                    -1.10488    0.19502
## PopulationStrawberry53                                    -0.61495    0.17821
## Test_environmentCherry                                     0.19537    0.20653
## Test_environmentStrawberry                                -0.46864    0.16128
## SA1                                                        0.33610    0.16217
## EggScore2                                                  0.71482    0.12567
## EggScore3                                                  1.00659    0.13193
## EggScore4                                                  1.12775    0.14440
## Test_environmentBlackberry:Original_environmentCherry      0.11735    0.21231
## Test_environmentCherry:Original_environmentCherry         -0.22493    0.25849
## Test_environmentStrawberry:Original_environmentCherry           NA         NA
## Test_environmentBlackberry:Original_environmentStrawberry  0.23700    0.28027
## Test_environmentCherry:Original_environmentStrawberry           NA         NA
## Test_environmentStrawberry:Original_environmentStrawberry       NA         NA
##                                                           t value Pr(>|t|)    
## (Intercept)                                                10.113  < 2e-16 ***
## PopulationBlackberry32                                    -12.591  < 2e-16 ***
## PopulationBlackberry33                                     -1.129 0.259492    
## PopulationBlackberry34                                     -4.678 3.62e-06 ***
## PopulationBlackberry35                                     -6.162 1.36e-09 ***
## PopulationBlackberry36                                     -1.872 0.061773 .  
## PopulationBlackberry37                                     -3.545 0.000425 ***
## PopulationBlackberry38                                     -4.818 1.86e-06 ***
## PopulationBlackberry39                                     -2.249 0.024885 *  
## PopulationBlackberry40                                     -0.464 0.642637    
## PopulationBlackberry43                                     -1.658 0.097822 .  
## PopulationBlackberry44                                     -0.166 0.867882    
## PopulationBlackberry45                                     -3.302 0.001021 ** 
## PopulationCherry103                                         0.557 0.577725    
## PopulationCherry104                                        -4.827 1.79e-06 ***
## PopulationCherry3                                         -10.533  < 2e-16 ***
## PopulationCherry47                                         -6.138 1.57e-09 ***
## PopulationCherry50                                          2.255 0.024507 *  
## PopulationCherry51                                          1.210 0.226659    
## PopulationCherry52                                          1.027 0.305093    
## PopulationCherry6                                          -7.188 2.08e-12 ***
## PopulationCherry7                                          -4.059 5.63e-05 ***
## PopulationStrawberry42                                     -0.072 0.942481    
## PopulationStrawberry44                                     -5.666 2.33e-08 ***
## PopulationStrawberry53                                     -3.451 0.000601 ***
## Test_environmentCherry                                      0.946 0.344561    
## Test_environmentStrawberry                                 -2.906 0.003806 ** 
## SA1                                                         2.072 0.038676 *  
## EggScore2                                                   5.688 2.06e-08 ***
## EggScore3                                                   7.630 9.98e-14 ***
## EggScore4                                                   7.810 2.78e-14 ***
## Test_environmentBlackberry:Original_environmentCherry       0.553 0.580680    
## Test_environmentCherry:Original_environmentCherry          -0.870 0.384579    
## Test_environmentStrawberry:Original_environmentCherry          NA       NA    
## Test_environmentBlackberry:Original_environmentStrawberry   0.846 0.398124    
## Test_environmentCherry:Original_environmentStrawberry          NA       NA    
## Test_environmentStrawberry:Original_environmentStrawberry      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6453 on 568 degrees of freedom
## Multiple R-squared:  0.6435, Adjusted R-squared:  0.6228 
## F-statistic: 31.07 on 33 and 568 DF,  p-value: < 2.2e-16
anova(m2)
## Analysis of Variance Table
## 
## Response: log(Nb_adults + 1)
##                                        Df Sum Sq Mean Sq F value    Pr(>F)    
## Population                             24 323.52  13.480 32.3685 < 2.2e-16 ***
## Test_environment                        2  69.66  34.828 83.6282 < 2.2e-16 ***
## SA                                      1   3.21   3.215  7.7195  0.005644 ** 
## EggScore                                3  30.19  10.063 24.1633 9.951e-15 ***
## Test_environment:Original_environment   3   0.36   0.121  0.2916  0.831493    
## Residuals                             568 236.55   0.416                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## F test for SA
(Fratio_int = (anova(m2)[3,2]/anova(m2)[5,2])/(1/anova(m2)[5,1]))
## [1] 26.47414
(pvalue_int = 1 - pf(Fratio_int, 1, anova(m2)[5,1])) #the correct test (see equation D7 in Appendix D of the paper)
## [1] 0.01422731
## Compute R2 = MS Interaction model without SA - MS Interaction model with SA / MS Interaction model without SA
(rsq2 <- 1-anova(m2)[5, 3]/((anova(m2)[3, 2]+anova(m2)[5, 2])/(anova(m2)[3, 1]+anova(m2)[5, 1])))
## [1] 0.8642878
## Compare with EggScoreFive
m3 <- lm(log(Nb_adults+1) ~ Population + Test_environment + SA + Original_environment:Test_environment+EggScoreFive, data=data[data$Generation=="G2",])

## Compare with EggScoreSmall
m4 <- lm(log(Nb_adults+1) ~ Population + Test_environment + SA + Original_environment:Test_environment+EggScoreSmall, data=data[data$Generation=="G2",])

## No egg information
m0 <- lm(log(Nb_adults+1) ~ Population + Test_environment + SA + Original_environment:Test_environment, data=data[data$Generation=="G2",])

model.sel(m0, m1, m2, m3, m4)
## Model selection table 
##     (Int) Ppl SA Tst_env Org_env:Tst_env log(Nb_egg+1) EgS ESF ESS df   logLik
## m4 2.3450   +  +       +               +                         + 37 -567.123
## m2 2.4300   +  +       +               +                 +         35 -573.037
## m3 2.4380   +  +       +               +                     +     36 -572.526
## m1 0.8623   +  +       +               +        0.5107             33 -576.218
## m0 3.2020   +  +       +               +                           32 -608.489
##      AICc delta weight
## m4 1213.2  0.00  0.952
## m2 1220.5  7.29  0.025
## m3 1221.8  8.54  0.013
## m1 1222.4  9.15  0.010
## m0 1284.7 71.46  0.000
## Models ranked by AICc(x)

2.4.3 Power Analysis

m2 <- lm(log(Nb_adults+1) ~ -1 + Population + Test_environment + Original_environment:Test_environment+EggScore, data=data[data$Generation=="G2",])
coef(m2)

simdata <- function(seed=1, nrep = 2, npop = 3, SAincrease=0){
  set.seed(seed)
  ## Sample populations
  PopBlackberry <- sort(sample(levels(data$Population)[grep("Blackberry", levels(data$Population))], size=npop, replace = TRUE))
  PopCherry <- sort(sample(levels(data$Population)[grep("Cherry", levels(data$Population))], size=npop, replace = TRUE))
  PopStrawberry <- sort(sample(levels(data$Population)[grep("Strawberry", levels(data$Population))], size=npop, replace = TRUE))
  
  ## Create dataset
  newdata <- expand.grid(Pop_new_name=1:(3*npop), Test_environment=levels(data$Test_environment), EggScore=factor(1:4), Rep=1:nrep)
  newdata$Original_environment <- rep(c("Blackberry", "Cherry", "Strawberry"), each=npop)
  
  ## Add original names to get the right coef in the model
  newdata$Pop_new_name <- as.factor(newdata$Pop_new_name)
  newdata$Population <- newdata$Pop_new_name
  levels(newdata$Population) <- c(PopBlackberry, PopCherry, PopStrawberry)
  
  ## Get design matrix
  mat <- model.matrix(~-1 + Population + Test_environment + Original_environment:Test_environment+EggScore, data=newdata)
  
  ## Check that the order of the coefficients is the same
  data.frame(colnames(mat), names(coef(m2))[names(coef(m2))%in%colnames(mat)])
  
  coefsim <- coef(m2)[names(coef(m2))%in%colnames(mat)]
  coefsim <- ifelse(is.na(coefsim), 0, coefsim)
  
  ## Increase fitness in sympatric environment
  coefsim[grepl("PopulationBlackberry|Test_environmentCherry:Original_environmentCherry|Test_environmentStrawberry:Original_environmentStrawberry", names(coefsim))] <- coefsim[grepl("PopulationBlackberry|Test_environmentCherry:Original_environmentCherry|Test_environmentStrawberry:Original_environmentStrawberry", names(coefsim))]+SAincrease
  
  newdata$resp <- mat%*%coefsim + rnorm(n=nrow(newdata), mean = 0, sd = sigma(m2))
  
  m3 <- lm(resp ~ -1+Population + Test_environment + Original_environment:Test_environment+EggScore, data=newdata)
  data.frame(coefsim, coef(m3))
  
  newdata$SA <- as.factor(ifelse(newdata$Original_environment==newdata$Test_environment, 1, 0))
  
  m4 <- lm(resp ~ Pop_new_name + Test_environment + SA + Original_environment:Test_environment+EggScore, data=newdata)
  summary(m4)
  anova(m4)
    
  ## F test for SA
  (Fratio_int = (anova(m4)[3,2]/anova(m4)[5,2])/(1/anova(m4)[5,1]))
  (pvalue_int = 1 - pf(Fratio_int, 1, anova(m4)[5,1])) #the correct test (see equation D7 in Appendix D of the paper)
  ## rsquare with single model
  (rsq2 <- 1-anova(m4)[5, 3]/((anova(m4)[3, 2]+anova(m4)[5, 2])/(anova(m4)[3, 1]+anova(m4)[5, 1])))
  #print(rsq2)
  return(pvalue_int)
}

simdata(seed=2, nrep = 30, npop = 3)

simpval <- sapply(1:500, simdata, nrep = 10, npop = 3)
mean(simpval)
## Compute power
mean(simpval<0.05)
hist(simpval)

simpval10rep <- sapply(1:500, simdata, nrep = 10, npop = 3)
simpval15rep <- sapply(1:500, simdata, nrep = 15, npop = 3)
simpval20rep <- sapply(1:500, simdata, nrep = 20, npop = 3)
simpval25rep <- sapply(1:500, simdata, nrep = 25, npop = 3)
simpval30rep <- sapply(1:500, simdata, nrep = 30, npop = 3)
simpval35rep <- sapply(1:500, simdata, nrep = 35, npop = 3)
simpval40rep <- sapply(1:500, simdata, nrep = 40, npop = 3)
simpval50rep <- sapply(1:500, simdata, nrep = 50, npop = 3)
simpval60rep <- sapply(1:500, simdata, nrep = 60, npop = 3)
simpval80rep <- sapply(1:500, simdata, nrep = 80, npop = 3)
simpval100rep <- sapply(1:500, simdata, nrep = 100, npop = 3)

## Compute power
val <- data.frame(simpval10rep, simpval15rep, simpval20rep, simpval25rep, simpval30rep, simpval35rep, simpval40rep, simpval50rep, simpval60rep, simpval80rep, simpval100rep)

computepower <- function(x){
  return(mean(x<0.05))
}

statpower <- apply(val, 2, computepower)

plot(c(10, 15, 20, 25, 30, 35, 40, 50, 60, 80, 100), statpower, las=1, xlab="Number of replicates", ylab="Statistical Power")

2.5 Emergence rate

2.5.1 Explore data

library(lattice)
xyplot(Emergence_Rate~Nb_eggs|Original_environment*Test_environment, data=data)

xyplot(Emergence_Rate~EggScore|Original_environment*Test_environment, data=data)

bwplot(Emergence_Rate~EggScore|Original_environment*Test_environment, data=data)

bwplot(Emergence_Rate~EggScoreFive|Original_environment*Test_environment, data=data)

bwplot(Emergence_Rate~EggScoreSmall|Original_environment*Test_environment, data=data)

AvEmergenceRate <- tapply(data$Emergence_Rate, list(data$EggScoreFive, data$Original_environment, data$Test_environment), mean)
tapply(data$Emergence_Rate, list(data$EggScoreFive, data$Original_environment, data$Test_environment), length)
## , , Blackberry
## 
##   Blackberry Cherry Strawberry
## 1         11      1          3
## 2         43     14         13
## 3         40     17         20
## 4         10     10          8
## 5          3      4          2
## 
## , , Cherry
## 
##   Blackberry Cherry Strawberry
## 1          1      2         NA
## 2         30     10         17
## 3         50     20         21
## 4         23     11          9
## 5          2      7         NA
## 
## , , Strawberry
## 
##   Blackberry Cherry Strawberry
## 1          7      3          6
## 2         55     22         13
## 3         38     16         16
## 4          7      6          9
## 5         NA     NA          2
AvEmergenceRate[, , "Blackberry"][, "Cherry"]/AvEmergenceRate[, , "Blackberry"][, "Blackberry"]
##         1         2         3         4         5 
## 0.2722049 0.8010658 0.7437340 0.9735309 0.1709919
AvEmergenceRate[, , "Cherry"][, "Blackberry"]/AvEmergenceRate[, , "Cherry"][, "Cherry"]
##         1         2         3         4         5 
## 0.0000000 0.9000996 1.1452193 0.8514873 1.4918117
AvEmergenceRate[, , "Blackberry"][, "Cherry"]/AvEmergenceRate[, , "Cherry"][, "Cherry"]
##         1         2         3         4         5 
## 0.8525097 0.9696583 1.1359488 0.9037047 0.4651315

2.5.2 ANR analysis

dataG2 <- data[data$Generation=="G2",]

m0 <- lm(asin(sqrt(Emergence_Rate)) ~ Population + Test_environment, data=dataG2)

dataG2$res <- residuals(m0)

m0 <- lm(res ~ -1 + Original_environment : Test_environment, data=dataG2)

summary(m0)
## 
## Call:
## lm(formula = res ~ -1 + Original_environment:Test_environment, 
##     data = dataG2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.41850 -0.09344 -0.00716  0.08081  0.93906 
## 
## Coefficients:
##                                                            Estimate Std. Error
## Original_environmentBlackberry:Test_environmentBlackberry  0.025940   0.014200
## Original_environmentCherry:Test_environmentBlackberry     -0.033385   0.021657
## Original_environmentStrawberry:Test_environmentBlackberry -0.026953   0.021657
## Original_environmentBlackberry:Test_environmentCherry     -0.009440   0.014267
## Original_environmentCherry:Test_environmentCherry          0.029225   0.020773
## Original_environmentStrawberry:Test_environmentCherry     -0.009799   0.021426
## Original_environmentBlackberry:Test_environmentStrawberry -0.016587   0.014200
## Original_environmentCherry:Test_environmentStrawberry      0.001584   0.021426
## Original_environmentStrawberry:Test_environmentStrawberry  0.036965   0.021657
##                                                           t value Pr(>|t|)  
## Original_environmentBlackberry:Test_environmentBlackberry   1.827   0.0682 .
## Original_environmentCherry:Test_environmentBlackberry      -1.541   0.1237  
## Original_environmentStrawberry:Test_environmentBlackberry  -1.245   0.2138  
## Original_environmentBlackberry:Test_environmentCherry      -0.662   0.5084  
## Original_environmentCherry:Test_environmentCherry           1.407   0.1600  
## Original_environmentStrawberry:Test_environmentCherry      -0.457   0.6476  
## Original_environmentBlackberry:Test_environmentStrawberry  -1.168   0.2432  
## Original_environmentCherry:Test_environmentStrawberry       0.074   0.9411  
## Original_environmentStrawberry:Test_environmentStrawberry   1.707   0.0884 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1469 on 593 degrees of freedom
## Multiple R-squared:  0.02334,    Adjusted R-squared:  0.008517 
## F-statistic: 1.575 on 9 and 593 DF,  p-value: 0.1193
m0 <- lm(asin(sqrt(Emergence_Rate)) ~ -1+ Population + Test_environment + Original_environment : Test_environment, data=dataG2)
summary(m0)
## 
## Call:
## lm(formula = asin(sqrt(Emergence_Rate)) ~ -1 + Population + Test_environment + 
##     Original_environment:Test_environment, data = dataG2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.41945 -0.09246 -0.00672  0.08070  0.93802 
## 
## Coefficients: (2 not defined because of singularities)
##                                                           Estimate Std. Error
## PopulationBlackberry31                                     0.63332    0.03114
## PopulationBlackberry32                                     0.27700    0.02673
## PopulationBlackberry33                                     0.57728    0.02542
## PopulationBlackberry34                                     0.41343    0.04653
## PopulationBlackberry35                                     0.43325    0.02511
## PopulationBlackberry36                                     0.52676    0.03841
## PopulationBlackberry37                                     0.52330    0.02862
## PopulationBlackberry38                                     0.31701    0.06781
## PopulationBlackberry39                                     0.51650    0.03277
## PopulationBlackberry40                                     0.66341    0.02978
## PopulationBlackberry43                                     0.63278    0.04895
## PopulationBlackberry44                                     0.60013    0.02978
## PopulationBlackberry45                                     0.40441    0.07536
## PopulationCherry103                                        0.69122    0.04398
## PopulationCherry104                                        0.42678    0.04646
## PopulationCherry3                                          0.16504    0.05681
## PopulationCherry47                                         0.38239    0.03622
## PopulationCherry50                                         0.77006    0.05241
## PopulationCherry51                                         0.92529    0.15428
## PopulationCherry52                                         0.78469    0.08027
## PopulationCherry6                                          0.24721    0.06208
## PopulationCherry7                                          0.45095    0.03677
## PopulationStrawberry42                                     0.67561    0.03518
## PopulationStrawberry44                                     0.42915    0.03860
## PopulationStrawberry53                                     0.57813    0.03315
## Test_environmentCherry                                    -0.07930    0.02054
## Test_environmentStrawberry                                -0.18264    0.02049
## Test_environmentBlackberry:Original_environmentCherry     -0.07859    0.03741
## Test_environmentCherry:Original_environmentCherry          0.02052    0.03673
## Test_environmentStrawberry:Original_environmentCherry           NA         NA
## Test_environmentBlackberry:Original_environmentStrawberry -0.10655    0.03733
## Test_environmentCherry:Original_environmentStrawberry     -0.05393    0.03722
## Test_environmentStrawberry:Original_environmentStrawberry       NA         NA
##                                                           t value Pr(>|t|)    
## PopulationBlackberry31                                     20.335  < 2e-16 ***
## PopulationBlackberry32                                     10.362  < 2e-16 ***
## PopulationBlackberry33                                     22.709  < 2e-16 ***
## PopulationBlackberry34                                      8.885  < 2e-16 ***
## PopulationBlackberry35                                     17.257  < 2e-16 ***
## PopulationBlackberry36                                     13.715  < 2e-16 ***
## PopulationBlackberry37                                     18.285  < 2e-16 ***
## PopulationBlackberry38                                      4.675 3.68e-06 ***
## PopulationBlackberry39                                     15.763  < 2e-16 ***
## PopulationBlackberry40                                     22.276  < 2e-16 ***
## PopulationBlackberry43                                     12.927  < 2e-16 ***
## PopulationBlackberry44                                     20.151  < 2e-16 ***
## PopulationBlackberry45                                      5.366 1.17e-07 ***
## PopulationCherry103                                        15.718  < 2e-16 ***
## PopulationCherry104                                         9.187  < 2e-16 ***
## PopulationCherry3                                           2.905 0.003814 ** 
## PopulationCherry47                                         10.556  < 2e-16 ***
## PopulationCherry50                                         14.693  < 2e-16 ***
## PopulationCherry51                                          5.997 3.56e-09 ***
## PopulationCherry52                                          9.776  < 2e-16 ***
## PopulationCherry6                                           3.982 7.72e-05 ***
## PopulationCherry7                                          12.265  < 2e-16 ***
## PopulationStrawberry42                                     19.205  < 2e-16 ***
## PopulationStrawberry44                                     11.117  < 2e-16 ***
## PopulationStrawberry53                                     17.439  < 2e-16 ***
## Test_environmentCherry                                     -3.861 0.000126 ***
## Test_environmentStrawberry                                 -8.914  < 2e-16 ***
## Test_environmentBlackberry:Original_environmentCherry      -2.101 0.036095 *  
## Test_environmentCherry:Original_environmentCherry           0.559 0.576564    
## Test_environmentStrawberry:Original_environmentCherry          NA       NA    
## Test_environmentBlackberry:Original_environmentStrawberry  -2.854 0.004474 ** 
## Test_environmentCherry:Original_environmentStrawberry      -1.449 0.147919    
## Test_environmentStrawberry:Original_environmentStrawberry      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1497 on 571 degrees of freedom
## Multiple R-squared:  0.9001, Adjusted R-squared:  0.8947 
## F-statistic:   166 on 31 and 571 DF,  p-value: < 2.2e-16
emmeans::emmeans(m0, list(pairwise ~ Original_environment:Test_environment), adjust = "tukey") #
## Warning in (function (object, at, cov.reduce = mean, cov.keep = get_emm_option("cov.keep"), : There are unevaluated constants in the response formula
## Auto-detection of the response transformation may be incorrect
## NOTE: A nesting structure was detected in the fitted model:
##     Population %in% (Test_environment*Original_environment), Original_environment %in% Test_environment
## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates
## $`emmeans of Original_environment, Test_environment`
##  Original_environment Test_environment emmean     SE  df lower.CL upper.CL
##  Blackberry           Blackberry        0.501 0.0161 571    0.470    0.533
##  Cherry               Blackberry        0.460 0.0273 571    0.406    0.513
##  Strawberry           Blackberry        0.454 0.0224 571    0.410    0.498
##  Blackberry           Cherry            0.422 0.0163 571    0.390    0.454
##  Cherry               Cherry            0.431 0.0239 571    0.384    0.478
##  Strawberry           Cherry            0.428 0.0222 571    0.384    0.471
##  Blackberry           Strawberry        0.319 0.0165 571    0.286    0.351
##  Cherry               Strawberry        0.307 0.0242 571    0.260    0.355
##  Strawberry           Strawberry        0.378 0.0224 571    0.334    0.422
## 
## Results are averaged over the levels of: Population 
## Results are given on the asin(sqrt(mu)) (not the response) scale. 
## Confidence level used: 0.95 
## 
## $`pairwise differences of Original_environment, Test_environment`
##  contrast                                      estimate     SE  df t.ratio
##  Blackberry Blackberry - Cherry Blackberry      0.04184 0.0317 571  1.319 
##  Blackberry Blackberry - Strawberry Blackberry  0.04702 0.0276 571  1.704 
##  Blackberry Blackberry - Blackberry Cherry      0.07930 0.0205 571  3.861 
##  Blackberry Blackberry - Cherry Cherry          0.07042 0.0289 571  2.441 
##  Blackberry Blackberry - Strawberry Cherry      0.07370 0.0274 571  2.688 
##  Blackberry Blackberry - Blackberry Strawberry  0.18264 0.0205 571  8.914 
##  Blackberry Blackberry - Cherry Strawberry      0.19428 0.0291 571  6.683 
##  Blackberry Blackberry - Strawberry Strawberry  0.12311 0.0276 571  4.460 
##  Cherry Blackberry - Strawberry Blackberry      0.00518 0.0353 571  0.147 
##  Cherry Blackberry - Blackberry Cherry          0.03746 0.0318 571  1.177 
##  Cherry Blackberry - Cherry Cherry              0.02858 0.0336 571  0.852 
##  Cherry Blackberry - Strawberry Cherry          0.03185 0.0352 571  0.905 
##  Cherry Blackberry - Blackberry Strawberry      0.14080 0.0319 571  4.414 
##  Cherry Blackberry - Cherry Strawberry          0.15244 0.0340 571  4.488 
##  Cherry Blackberry - Strawberry Strawberry      0.08127 0.0353 571  2.300 
##  Strawberry Blackberry - Blackberry Cherry      0.03228 0.0277 571  1.164 
##  Strawberry Blackberry - Cherry Cherry          0.02340 0.0328 571  0.713 
##  Strawberry Blackberry - Strawberry Cherry      0.02667 0.0310 571  0.859 
##  Strawberry Blackberry - Blackberry Strawberry  0.13562 0.0278 571  4.875 
##  Strawberry Blackberry - Cherry Strawberry      0.14726 0.0330 571  4.464 
##  Strawberry Blackberry - Strawberry Strawberry  0.07609 0.0312 571  2.438 
##  Blackberry Cherry - Cherry Cherry             -0.00888 0.0290 571 -0.306 
##  Blackberry Cherry - Strawberry Cherry         -0.00560 0.0275 571 -0.203 
##  Blackberry Cherry - Blackberry Strawberry      0.10334 0.0205 571  5.033 
##  Blackberry Cherry - Cherry Strawberry          0.11498 0.0292 571  3.938 
##  Blackberry Cherry - Strawberry Strawberry      0.04381 0.0277 571  1.580 
##  Cherry Cherry - Strawberry Cherry              0.00328 0.0326 571  0.100 
##  Cherry Cherry - Blackberry Strawberry          0.11222 0.0291 571  3.862 
##  Cherry Cherry - Cherry Strawberry              0.12386 0.0304 571  4.068 
##  Cherry Cherry - Strawberry Strawberry          0.05269 0.0328 571  1.607 
##  Strawberry Cherry - Blackberry Strawberry      0.10895 0.0276 571  3.943 
##  Strawberry Cherry - Cherry Strawberry          0.12059 0.0328 571  3.673 
##  Strawberry Cherry - Strawberry Strawberry      0.04941 0.0310 571  1.592 
##  Blackberry Strawberry - Cherry Strawberry      0.01164 0.0293 571  0.398 
##  Blackberry Strawberry - Strawberry Strawberry -0.05953 0.0278 571 -2.140 
##  Cherry Strawberry - Strawberry Strawberry     -0.07117 0.0330 571 -2.158 
##  p.value
##  0.9254 
##  0.7443 
##  0.0040 
##  0.2640 
##  0.1543 
##  <.0001 
##  <.0001 
##  0.0003 
##  1.0000 
##  0.9611 
##  0.9952 
##  0.9927 
##  0.0004 
##  0.0003 
##  0.3440 
##  0.9636 
##  0.9986 
##  0.9948 
##  <.0001 
##  0.0003 
##  0.2655 
##  1.0000 
##  1.0000 
##  <.0001 
##  0.0030 
##  0.8157 
##  1.0000 
##  0.0040 
##  0.0018 
##  0.8011 
##  0.0029 
##  0.0080 
##  0.8093 
##  1.0000 
##  0.4469 
##  0.4352 
## 
## Results are averaged over the levels of: Population 
## Note: contrasts are still on the asin.sqrt scale 
## P value adjustment: tukey method for comparing a family of 9 estimates
dataG2$Nb_adults_corr <- dataG2$Nb_adults
dataG2$Nb_adults_corr[dataG2$Nb_adults>dataG2$Nb_eggs] <- dataG2$Nb_eggs[dataG2$Nb_adults>dataG2$Nb_eggs] 

m0 <- glm(cbind(Nb_adults_corr, Nb_eggs-Nb_adults_corr) ~ -1+ Population + Test_environment, family=binomial, data=dataG2)

hist(residuals(m0))

m0 <- lm(residuals(m0) ~ -1 + Original_environment : Test_environment, data=dataG2)
summary(m0)
## 
## Call:
## lm(formula = residuals(m0) ~ -1 + Original_environment:Test_environment, 
##     data = dataG2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.3245 -1.8429 -0.1601  1.6310 12.7946 
## 
## Coefficients:
##                                                           Estimate Std. Error
## Original_environmentBlackberry:Test_environmentBlackberry   0.2566     0.2570
## Original_environmentCherry:Test_environmentBlackberry      -0.7885     0.3920
## Original_environmentStrawberry:Test_environmentBlackberry  -0.4115     0.3920
## Original_environmentBlackberry:Test_environmentCherry      -0.2386     0.2582
## Original_environmentCherry:Test_environmentCherry           0.5971     0.3760
## Original_environmentStrawberry:Test_environmentCherry       0.0459     0.3878
## Original_environmentBlackberry:Test_environmentStrawberry  -0.5650     0.2570
## Original_environmentCherry:Test_environmentStrawberry      -0.2881     0.3878
## Original_environmentStrawberry:Test_environmentStrawberry   0.4425     0.3920
##                                                           t value Pr(>|t|)  
## Original_environmentBlackberry:Test_environmentBlackberry   0.998   0.3186  
## Original_environmentCherry:Test_environmentBlackberry      -2.011   0.0447 *
## Original_environmentStrawberry:Test_environmentBlackberry  -1.050   0.2943  
## Original_environmentBlackberry:Test_environmentCherry      -0.924   0.3559  
## Original_environmentCherry:Test_environmentCherry           1.588   0.1129  
## Original_environmentStrawberry:Test_environmentCherry       0.118   0.9058  
## Original_environmentBlackberry:Test_environmentStrawberry  -2.198   0.0283 *
## Original_environmentCherry:Test_environmentStrawberry      -0.743   0.4579  
## Original_environmentStrawberry:Test_environmentStrawberry   1.129   0.2594  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.659 on 593 degrees of freedom
## Multiple R-squared:  0.02658,    Adjusted R-squared:  0.0118 
## F-statistic: 1.799 on 9 and 593 DF,  p-value: 0.06551
emmeans::emmeans(m0, list(pairwise ~ Original_environment:Test_environment), adjust = "tukey") #
## NOTE: A nesting structure was detected in the fitted model:
##     Original_environment %in% (Original_environment*Test_environment), Test_environment %in% (Original_environment*Test_environment)
## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates
## $`emmeans of Original_environment, Test_environment`
##  Original_environment Test_environment  emmean    SE  df lower.CL upper.CL
##  Blackberry           Blackberry        0.2566 0.257 593   -0.248   0.7614
##  Cherry               Blackberry       -0.7885 0.392 593   -1.558  -0.0185
##  Strawberry           Blackberry       -0.4115 0.392 593   -1.181   0.3585
##  Blackberry           Cherry           -0.2386 0.258 593   -0.746   0.2686
##  Cherry               Cherry            0.5971 0.376 593   -0.141   1.3355
##  Strawberry           Cherry            0.0459 0.388 593   -0.716   0.8076
##  Blackberry           Strawberry       -0.5650 0.257 593   -1.070  -0.0602
##  Cherry               Strawberry       -0.2881 0.388 593   -1.050   0.4736
##  Strawberry           Strawberry        0.4425 0.392 593   -0.327   1.2125
## 
## Results are given on the residuals (not the response) scale. 
## Confidence level used: 0.95 
## 
## $`pairwise differences of Original_environment, Test_environment`
##  contrast                                      estimate    SE  df t.ratio
##  Blackberry Blackberry - Cherry Blackberry       1.0450 0.469 593  2.229 
##  Blackberry Blackberry - Strawberry Blackberry   0.6680 0.469 593  1.425 
##  Blackberry Blackberry - Blackberry Cherry       0.4952 0.364 593  1.359 
##  Blackberry Blackberry - Cherry Cherry          -0.3405 0.455 593 -0.748 
##  Blackberry Blackberry - Strawberry Cherry       0.2107 0.465 593  0.453 
##  Blackberry Blackberry - Blackberry Strawberry   0.8216 0.364 593  2.260 
##  Blackberry Blackberry - Cherry Strawberry       0.5447 0.465 593  1.171 
##  Blackberry Blackberry - Strawberry Strawberry  -0.1860 0.469 593 -0.397 
##  Cherry Blackberry - Strawberry Blackberry      -0.3770 0.554 593 -0.680 
##  Cherry Blackberry - Blackberry Cherry          -0.5499 0.469 593 -1.171 
##  Cherry Blackberry - Cherry Cherry              -1.3855 0.543 593 -2.551 
##  Cherry Blackberry - Strawberry Cherry          -0.8344 0.551 593 -1.513 
##  Cherry Blackberry - Blackberry Strawberry      -0.2235 0.469 593 -0.477 
##  Cherry Blackberry - Cherry Strawberry          -0.5004 0.551 593 -0.907 
##  Cherry Blackberry - Strawberry Strawberry      -1.2310 0.554 593 -2.220 
##  Strawberry Blackberry - Blackberry Cherry      -0.1729 0.469 593 -0.368 
##  Strawberry Blackberry - Cherry Cherry          -1.0085 0.543 593 -1.857 
##  Strawberry Blackberry - Strawberry Cherry      -0.4574 0.551 593 -0.829 
##  Strawberry Blackberry - Blackberry Strawberry   0.1535 0.469 593  0.327 
##  Strawberry Blackberry - Cherry Strawberry      -0.1234 0.551 593 -0.224 
##  Strawberry Blackberry - Strawberry Strawberry  -0.8540 0.554 593 -1.540 
##  Blackberry Cherry - Cherry Cherry              -0.8356 0.456 593 -1.832 
##  Blackberry Cherry - Strawberry Cherry          -0.2845 0.466 593 -0.611 
##  Blackberry Cherry - Blackberry Strawberry       0.3264 0.364 593  0.896 
##  Blackberry Cherry - Cherry Strawberry           0.0495 0.466 593  0.106 
##  Blackberry Cherry - Strawberry Strawberry      -0.6811 0.469 593 -1.451 
##  Cherry Cherry - Strawberry Cherry               0.5512 0.540 593  1.020 
##  Cherry Cherry - Blackberry Strawberry           1.1620 0.455 593  2.551 
##  Cherry Cherry - Cherry Strawberry               0.8851 0.540 593  1.639 
##  Cherry Cherry - Strawberry Strawberry           0.1545 0.543 593  0.284 
##  Strawberry Cherry - Blackberry Strawberry       0.6109 0.465 593  1.313 
##  Strawberry Cherry - Cherry Strawberry           0.3340 0.548 593  0.609 
##  Strawberry Cherry - Strawberry Strawberry      -0.3966 0.551 593 -0.719 
##  Blackberry Strawberry - Cherry Strawberry      -0.2769 0.465 593 -0.595 
##  Blackberry Strawberry - Strawberry Strawberry  -1.0075 0.469 593 -2.149 
##  Cherry Strawberry - Strawberry Strawberry      -0.7306 0.551 593 -1.325 
##  p.value
##  0.3879 
##  0.8880 
##  0.9125 
##  0.9981 
##  1.0000 
##  0.3684 
##  0.9623 
##  1.0000 
##  0.9990 
##  0.9622 
##  0.2102 
##  0.8493 
##  0.9999 
##  0.9925 
##  0.3937 
##  1.0000 
##  0.6440 
##  0.9960 
##  1.0000 
##  1.0000 
##  0.8360 
##  0.6608 
##  0.9996 
##  0.9932 
##  1.0000 
##  0.8773 
##  0.9839 
##  0.2099 
##  0.7831 
##  1.0000 
##  0.9274 
##  0.9996 
##  0.9985 
##  0.9996 
##  0.4407 
##  0.9237 
## 
## Note: contrasts are still on the residuals scale 
## P value adjustment: tukey method for comparing a family of 9 estimates

2.5.3 Explore data with blue

m0 <- lm(Emergence_Rate~Original_environment*Test_environment, data=data[data$Generation=="G2",])

summary(m0)
## 
## Call:
## lm(formula = Emergence_Rate ~ Original_environment * Test_environment, 
##     data = data[data$Generation == "G2", ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.26484 -0.09978 -0.02467  0.07237  0.73516 
## 
## Coefficients:
##                                                           Estimate Std. Error
## (Intercept)                                                0.26484    0.01330
## Original_environmentCherry                                -0.07926    0.02426
## Original_environmentStrawberry                            -0.04419    0.02426
## Test_environmentCherry                                    -0.07253    0.01886
## Test_environmentStrawberry                                -0.13033    0.01881
## Original_environmentCherry:Test_environmentCherry          0.06711    0.03385
## Original_environmentStrawberry:Test_environmentCherry      0.04696    0.03421
## Original_environmentCherry:Test_environmentStrawberry      0.05509    0.03418
## Original_environmentStrawberry:Test_environmentStrawberry  0.07162    0.03431
##                                                           t value Pr(>|t|)    
## (Intercept)                                                19.909  < 2e-16 ***
## Original_environmentCherry                                 -3.267 0.001150 ** 
## Original_environmentStrawberry                             -1.822 0.069020 .  
## Test_environmentCherry                                     -3.846 0.000133 ***
## Test_environmentStrawberry                                 -6.928 1.12e-11 ***
## Original_environmentCherry:Test_environmentCherry           1.982 0.047891 *  
## Original_environmentStrawberry:Test_environmentCherry       1.373 0.170346    
## Original_environmentCherry:Test_environmentStrawberry       1.612 0.107534    
## Original_environmentStrawberry:Test_environmentStrawberry   2.088 0.037265 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1376 on 593 degrees of freedom
## Multiple R-squared:  0.1044, Adjusted R-squared:  0.09231 
## F-statistic:  8.64 on 8 and 593 DF,  p-value: 3.499e-11
tapply(data$Emergence_Rate[data$Generation=="G2"], list(data$Original_environment[data$Generation=="G2"], data$Test_environment[data$Generation=="G2"]), mean)
##            Blackberry    Cherry Strawberry
## Blackberry  0.2648444 0.1923140  0.1345134
## Cherry      0.1855892 0.1801664  0.1103531
## Strawberry  0.2206513 0.1950778  0.1619433
tapply(data$Emergence_Rate[data$Generation=="G2"], list(data$Original_environment[data$Generation=="G2"], data$Test_environment[data$Generation=="G2"]), var)
##            Blackberry      Cherry Strawberry
## Blackberry 0.03491514 0.013749525 0.01436498
## Cherry     0.02922031 0.019467766 0.01399160
## Strawberry 0.01555929 0.009464032 0.01139614
range(data$Nb_adults)
## [1]   0 108
## Check number of eggs and adults
tapply(data$Nb_eggs[data$Generation=="G2"], list(data$Population[data$Generation=="G2"], data$Test_environment[data$Generation=="G2"]), mean)
##              Blackberry    Cherry Strawberry
## Blackberry31  106.44444 130.33333  118.88889
## Blackberry32   61.23077  84.46154   74.00000
## Blackberry33  125.13333 141.28571   95.66667
## Blackberry34  105.75000  95.66667   83.50000
## Blackberry35   93.60000 123.73333   96.31250
## Blackberry36  118.00000 144.00000   92.83333
## Blackberry37  108.09091 114.45455  101.27273
## Blackberry38  125.50000 141.00000  143.00000
## Blackberry39  133.75000 146.00000  117.12500
## Blackberry40   88.30000  92.10000   80.40000
## Blackberry43   69.33333  90.00000   86.33333
## Blackberry44  129.20000 155.70000  129.70000
## Blackberry45   92.00000 113.00000  110.00000
## Cherry103     129.66667 118.00000   89.83333
## Cherry104      59.80000  61.83333   55.20000
## Cherry3       109.66667 143.66667   78.00000
## Cherry47      142.00000 162.71429  128.23077
## Cherry50      148.33333 185.00000  126.00000
## Cherry51      139.00000        NA         NA
## Cherry52      129.50000  96.00000  103.00000
## Cherry6       226.50000 203.50000  124.00000
## Cherry7       137.75000 118.69231   97.50000
## Strawberry42  140.46667 131.81250  146.66667
## Strawberry44  150.60000 142.30000  139.70000
## Strawberry53   89.00000  99.23810   73.09524
tapply(data$Nb_adults[data$Generation=="G2"], list(data$Population[data$Generation=="G2"], data$Test_environment[data$Generation=="G2"]), mean)
##              Blackberry    Cherry Strawberry
## Blackberry31  38.444444 36.777778 19.2222222
## Blackberry32   1.923077  6.307692  0.6923077
## Blackberry33  30.800000 34.000000 16.6666667
## Blackberry34  13.500000 17.666667  8.0000000
## Blackberry35  14.533333 23.533333  5.2500000
## Blackberry36  28.800000 23.166667 14.0000000
## Blackberry37  34.727273 24.090909  8.2727273
## Blackberry38  15.500000  7.000000  1.0000000
## Blackberry39  28.375000 26.750000 19.0000000
## Blackberry40  38.700000 17.300000 23.9000000
## Blackberry43  45.000000 13.750000 11.6666667
## Blackberry44  48.900000 27.300000 25.6000000
## Blackberry45  11.500000 30.000000  2.0000000
## Cherry103     35.666667 41.857143 20.8333333
## Cherry104      8.200000 10.833333  3.0000000
## Cherry3        1.000000  3.000000  0.0000000
## Cherry47      15.833333 18.571429  5.9230769
## Cherry50      63.000000 67.750000 46.0000000
## Cherry51      78.000000        NA         NA
## Cherry52      65.500000 39.000000 18.0000000
## Cherry6        5.500000 10.000000  1.3333333
## Cherry7       16.416667 16.076923  9.9166667
## Strawberry42  45.733333 35.000000 26.3333333
## Strawberry44  13.100000 13.400000  9.9000000
## Strawberry53  17.285714 17.000000 11.6190476
## Check for the presence of negative correlations
m1 <- lm(asin(sqrt(Emergence_Rate)) ~ Population + Test_environment, data=data[data$Generation=="G2",])
data$resid <- residuals(m1)

meanbypopbytestenv <- as.data.frame(tapply(data$resid[data$Generation=="G2"], list(data$Population[data$Generation=="G2"], data$Test_environment[data$Generation=="G2"]), mean))


## Cherry ~ Blackberry
plot(meanbypopbytestenv$Blackberry, meanbypopbytestenv$Cherry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), xlab="Blackberry", ylab="Cherry", pch=16)
legend("topright", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)

## Strawberry ~ Cherry
plot(meanbypopbytestenv$Cherry, meanbypopbytestenv$Strawberry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), xlab="Cherry", ylab="Strawberry", pch=16)
legend("bottomright", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)

## Strawberry ~ Blackberry
plot(meanbypopbytestenv$Blackberry~ meanbypopbytestenv$Strawberry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), xlab="Blackberry", ylab="Strawberry", pch=16)
legend("topright", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)

2.5.4 Explore data with BLUPS

m0 <- lmer(asin(sqrt(Emergence_Rate))~ Population + Test_environment + (0 + Test_environment | Population), data = data[data$Generation=="G2",])
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 2 negative eigenvalues
## Correlation between fixed effects
plot(coef(m1), fixef(m0))

summary(m0)
## Linear mixed model fit by REML ['lmerMod']
## Formula: asin(sqrt(Emergence_Rate)) ~ Population + Test_environment +  
##     (0 + Test_environment | Population)
##    Data: data[data$Generation == "G2", ]
## 
## REML criterion at convergence: -480.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4905 -0.6208 -0.0492  0.5313  5.5124 
## 
## Random effects:
##  Groups     Name                       Variance  Std.Dev. Corr     
##  Population Test_environmentBlackberry 0.0009946 0.03154           
##             Test_environmentCherry     0.0114799 0.10714  0.08     
##             Test_environmentStrawberry 0.0080802 0.08989  0.33 0.79
##  Residual                              0.0206939 0.14385           
## Number of obs: 602, groups:  Population, 25
## 
## Fixed effects:
##                            Estimate Std. Error t value
## (Intercept)                 0.65326    0.05211  12.536
## PopulationBlackberry32     -0.42343    0.07005  -6.045
## PopulationBlackberry33     -0.11104    0.06898  -1.610
## PopulationBlackberry34     -0.31902    0.08364  -3.814
## PopulationBlackberry35     -0.25882    0.06897  -3.753
## PopulationBlackberry36     -0.14211    0.07997  -1.777
## PopulationBlackberry37     -0.08778    0.07144  -1.229
## PopulationBlackberry38     -0.33401    0.09810  -3.405
## PopulationBlackberry39     -0.17511    0.07451  -2.350
## PopulationBlackberry40      0.03641    0.07230   0.504
## PopulationBlackberry43      0.11160    0.08780   1.271
## PopulationBlackberry44     -0.02326    0.07230  -0.322
## PopulationBlackberry45     -0.27941    0.10044  -2.782
## PopulationCherry103        -0.02717    0.07774  -0.350
## PopulationCherry104        -0.27661    0.08008  -3.454
## PopulationCherry3          -0.55496    0.08815  -6.296
## PopulationCherry47         -0.32855    0.07066  -4.649
## PopulationCherry50          0.05616    0.08738   0.643
## PopulationCherry51          0.19344    0.15622   1.238
## PopulationCherry52          0.09981    0.10044   0.994
## PopulationCherry6          -0.47735    0.09491  -5.030
## PopulationCherry7          -0.28451    0.07068  -4.025
## PopulationStrawberry42     -0.04642    0.06896  -0.673
## PopulationStrawberry44     -0.34414    0.07230  -4.760
## PopulationStrawberry53     -0.19139    0.06681  -2.865
## Test_environmentCherry     -0.04937    0.02799  -1.764
## Test_environmentStrawberry -0.14711    0.02380  -6.180
## 
## Correlation matrix not shown by default, as p = 27 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
## optimizer (nloptwrap) convergence code: 0 (OK)
## unable to evaluate scaled gradient
## Model failed to converge: degenerate  Hessian with 2 negative eigenvalues
BLUP <- ranef(m0)$Population

pairs(BLUP)

## Positive correlation betwen blups and blues, as expected
plot(BLUP$Test_environmentBlackberry, meanbypopbytestenv$Blackberry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), pch=16)
legend("topleft", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)

## Positive correlation betwen blups and blues, as expected
plot(BLUP$Test_environmentCherry, meanbypopbytestenv$Cherry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), pch=16)
legend("topleft", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)

## Positive correlation betwen blups and blues, as expected
plot(BLUP$Test_environmentStrawberry, meanbypopbytestenv$Strawberry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), pch=16)
legend("bottomright", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)

2.5.5 Analysis

data$SA <- as.factor(ifelse(data$Original_environment==data$Test_environment, 1, 0))

m1 <- lm(asin(sqrt(Emergence_Rate)) ~ Population + Test_environment + SA + Original_environment:Test_environment, data=data[data$Generation=="G2",])
summary(m1)
## 
## Call:
## lm(formula = asin(sqrt(Emergence_Rate)) ~ Population + Test_environment + 
##     SA + Original_environment:Test_environment, data = data[data$Generation == 
##     "G2", ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.41945 -0.09246 -0.00672  0.08070  0.93802 
## 
## Coefficients: (3 not defined because of singularities)
##                                                             Estimate Std. Error
## (Intercept)                                                0.5793941  0.0485461
## PopulationBlackberry32                                    -0.3563256  0.0374731
## PopulationBlackberry33                                    -0.0560407  0.0365927
## PopulationBlackberry34                                    -0.2198934  0.0535505
## PopulationBlackberry35                                    -0.2000730  0.0362889
## PopulationBlackberry36                                    -0.1065683  0.0463480
## PopulationBlackberry37                                    -0.1100221  0.0388418
## PopulationBlackberry38                                    -0.3163117  0.0729120
## PopulationBlackberry39                                    -0.1168236  0.0419913
## PopulationBlackberry40                                     0.0300903  0.0397061
## PopulationBlackberry43                                    -0.0005464  0.0554218
## PopulationBlackberry44                                    -0.0331930  0.0397061
## PopulationBlackberry45                                    -0.2289148  0.0802466
## PopulationCherry103                                        0.0578950  0.0498387
## PopulationCherry104                                       -0.2065488  0.0520413
## PopulationCherry3                                         -0.4682877  0.0614604
## PopulationCherry47                                        -0.2509357  0.0431538
## PopulationCherry50                                         0.1367322  0.0574191
## PopulationCherry51                                         0.2919698  0.1560561
## PopulationCherry52                                         0.1513627  0.0836261
## PopulationCherry6                                         -0.3861147  0.0663642
## PopulationCherry7                                         -0.1823742  0.0436110
## PopulationStrawberry42                                    -0.0116455  0.0419997
## PopulationStrawberry44                                    -0.2581011  0.0451386
## PopulationStrawberry53                                    -0.1091287  0.0405748
## Test_environmentCherry                                    -0.0253687  0.0472313
## Test_environmentStrawberry                                -0.1287121  0.0372242
## SA1                                                        0.0539301  0.0372219
## Test_environmentBlackberry:Original_environmentCherry     -0.0246612  0.0486347
## Test_environmentCherry:Original_environmentCherry         -0.0334106  0.0598113
## Test_environmentStrawberry:Original_environmentCherry             NA         NA
## Test_environmentBlackberry:Original_environmentStrawberry  0.0013057  0.0644075
## Test_environmentCherry:Original_environmentStrawberry             NA         NA
## Test_environmentStrawberry:Original_environmentStrawberry         NA         NA
##                                                           t value Pr(>|t|)    
## (Intercept)                                                11.935  < 2e-16 ***
## PopulationBlackberry32                                     -9.509  < 2e-16 ***
## PopulationBlackberry33                                     -1.531 0.126207    
## PopulationBlackberry34                                     -4.106 4.61e-05 ***
## PopulationBlackberry35                                     -5.513 5.34e-08 ***
## PopulationBlackberry36                                     -2.299 0.021847 *  
## PopulationBlackberry37                                     -2.833 0.004780 ** 
## PopulationBlackberry38                                     -4.338 1.70e-05 ***
## PopulationBlackberry39                                     -2.782 0.005580 ** 
## PopulationBlackberry40                                      0.758 0.448868    
## PopulationBlackberry43                                     -0.010 0.992137    
## PopulationBlackberry44                                     -0.836 0.403523    
## PopulationBlackberry45                                     -2.853 0.004493 ** 
## PopulationCherry103                                         1.162 0.245864    
## PopulationCherry104                                        -3.969 8.14e-05 ***
## PopulationCherry3                                          -7.619 1.07e-13 ***
## PopulationCherry47                                         -5.815 1.01e-08 ***
## PopulationCherry50                                          2.381 0.017578 *  
## PopulationCherry51                                          1.871 0.061866 .  
## PopulationCherry52                                          1.810 0.070822 .  
## PopulationCherry6                                          -5.818 9.92e-09 ***
## PopulationCherry7                                          -4.182 3.35e-05 ***
## PopulationStrawberry42                                     -0.277 0.781668    
## PopulationStrawberry44                                     -5.718 1.74e-08 ***
## PopulationStrawberry53                                     -2.690 0.007364 ** 
## Test_environmentCherry                                     -0.537 0.591397    
## Test_environmentStrawberry                                 -3.458 0.000585 ***
## SA1                                                         1.449 0.147919    
## Test_environmentBlackberry:Original_environmentCherry      -0.507 0.612301    
## Test_environmentCherry:Original_environmentCherry          -0.559 0.576653    
## Test_environmentStrawberry:Original_environmentCherry          NA       NA    
## Test_environmentBlackberry:Original_environmentStrawberry   0.020 0.983833    
## Test_environmentCherry:Original_environmentStrawberry          NA       NA    
## Test_environmentStrawberry:Original_environmentStrawberry      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1497 on 571 degrees of freedom
## Multiple R-squared:  0.493,  Adjusted R-squared:  0.4664 
## F-statistic: 18.51 on 30 and 571 DF,  p-value: < 2.2e-16
anova(m1)
## Analysis of Variance Table
## 
## Response: asin(sqrt(Emergence_Rate))
##                                        Df  Sum Sq Mean Sq F value    Pr(>F)    
## Population                             24 10.0931 0.42054 18.7710 < 2.2e-16 ***
## Test_environment                        2  2.0401 1.02004 45.5294 < 2.2e-16 ***
## SA                                      1  0.2890 0.28899 12.8989 0.0003571 ***
## Test_environment:Original_environment   3  0.0186 0.00621  0.2773 0.8417788    
## Residuals                             571 12.7926 0.02240                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## F test for SA
(Fratio_int = (anova(m1)[3,2]/anova(m1)[4,2])/(1/anova(m1)[4,1]))
## [1] 46.51257
(pvalue_int = 1 - pf(Fratio_int, 1, anova(m1)[4,1])) #the correct test (see equation D7 in Appendix D of the paper)
## [1] 0.006448859
## Compute Rsquare
(rsq2 <- 1-anova(m1)[4, 3]/((anova(m1)[3, 2]+anova(m1)[4, 2])/(anova(m1)[3, 1]+anova(m1)[4, 1])))
## [1] 0.9192124
## Correct for number of eggs
m2 <- lm(asin(sqrt(Emergence_Rate)) ~ Population + Test_environment + SA + Original_environment:Test_environment+log(Nb_eggs+1), data=data[data$Generation=="G2",])

## F test for SA
(Fratio_int = (anova(m2)[3,2]/anova(m2)[5,2])/(1/anova(m2)[5,1]))
## [1] 108.7824
(pvalue_int = 1 - pf(Fratio_int, 1, anova(m2)[5,1])) #the correct test (see equation D7 in Appendix D of the paper)
## [1] 0.001881239
## Compute Rsquare
(rsq2 <- 1-anova(m2)[5, 3]/((anova(m2)[3, 2]+anova(m2)[5, 2])/(anova(m2)[3, 1]+anova(m2)[5, 1])))
## [1] 0.9642162
m3 <- lm(asin(sqrt(Emergence_Rate)) ~ Population + Test_environment + SA + Original_environment:Test_environment+EggScore, data=data[data$Generation=="G2",])
summary(m3)
## 
## Call:
## lm(formula = asin(sqrt(Emergence_Rate)) ~ Population + Test_environment + 
##     SA + Original_environment:Test_environment + EggScore, data = data[data$Generation == 
##     "G2", ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.43011 -0.09039 -0.00357  0.07908  0.93919 
## 
## Coefficients: (3 not defined because of singularities)
##                                                            Estimate Std. Error
## (Intercept)                                                0.602001   0.055208
## PopulationBlackberry32                                    -0.377760   0.037911
## PopulationBlackberry33                                    -0.055998   0.036326
## PopulationBlackberry34                                    -0.229305   0.053200
## PopulationBlackberry35                                    -0.210254   0.036129
## PopulationBlackberry36                                    -0.107246   0.046039
## PopulationBlackberry37                                    -0.120148   0.038631
## PopulationBlackberry38                                    -0.321473   0.072756
## PopulationBlackberry39                                    -0.108157   0.041700
## PopulationBlackberry40                                     0.011493   0.039768
## PopulationBlackberry43                                    -0.023178   0.055325
## PopulationBlackberry44                                    -0.027083   0.039655
## PopulationBlackberry45                                    -0.237178   0.079741
## PopulationCherry103                                        0.045284   0.049575
## PopulationCherry104                                       -0.240476   0.052701
## PopulationCherry3                                         -0.480668   0.061042
## PopulationCherry47                                        -0.240715   0.042888
## PopulationCherry50                                         0.158899   0.057202
## PopulationCherry51                                         0.278850   0.154902
## PopulationCherry52                                         0.136201   0.083137
## PopulationCherry6                                         -0.367035   0.065953
## PopulationCherry7                                         -0.191145   0.043298
## PopulationStrawberry42                                    -0.008004   0.041722
## PopulationStrawberry44                                    -0.255111   0.044814
## PopulationStrawberry53                                    -0.136834   0.040950
## Test_environmentCherry                                     0.001684   0.047459
## Test_environmentStrawberry                                -0.117301   0.037060
## SA1                                                        0.069153   0.037266
## EggScore2                                                 -0.016370   0.028877
## EggScore3                                                 -0.036638   0.030316
## EggScore4                                                 -0.088070   0.033183
## Test_environmentBlackberry:Original_environmentCherry      0.002489   0.048787
## Test_environmentCherry:Original_environmentCherry         -0.044893   0.059398
## Test_environmentStrawberry:Original_environmentCherry            NA         NA
## Test_environmentBlackberry:Original_environmentStrawberry  0.029211   0.064404
## Test_environmentCherry:Original_environmentStrawberry            NA         NA
## Test_environmentStrawberry:Original_environmentStrawberry        NA         NA
##                                                           t value Pr(>|t|)    
## (Intercept)                                                10.904  < 2e-16 ***
## PopulationBlackberry32                                     -9.964  < 2e-16 ***
## PopulationBlackberry33                                     -1.542 0.123741    
## PopulationBlackberry34                                     -4.310 1.92e-05 ***
## PopulationBlackberry35                                     -5.820 9.86e-09 ***
## PopulationBlackberry36                                     -2.329 0.020184 *  
## PopulationBlackberry37                                     -3.110 0.001964 ** 
## PopulationBlackberry38                                     -4.419 1.19e-05 ***
## PopulationBlackberry39                                     -2.594 0.009741 ** 
## PopulationBlackberry40                                      0.289 0.772692    
## PopulationBlackberry43                                     -0.419 0.675418    
## PopulationBlackberry44                                     -0.683 0.494907    
## PopulationBlackberry45                                     -2.974 0.003061 ** 
## PopulationCherry103                                         0.913 0.361395    
## PopulationCherry104                                        -4.563 6.18e-06 ***
## PopulationCherry3                                          -7.874 1.75e-14 ***
## PopulationCherry47                                         -5.613 3.12e-08 ***
## PopulationCherry50                                          2.778 0.005653 ** 
## PopulationCherry51                                          1.800 0.072364 .  
## PopulationCherry52                                          1.638 0.101919    
## PopulationCherry6                                          -5.565 4.04e-08 ***
## PopulationCherry7                                          -4.415 1.21e-05 ***
## PopulationStrawberry42                                     -0.192 0.847934    
## PopulationStrawberry44                                     -5.693 2.01e-08 ***
## PopulationStrawberry53                                     -3.341 0.000888 ***
## Test_environmentCherry                                      0.035 0.971714    
## Test_environmentStrawberry                                 -3.165 0.001633 ** 
## SA1                                                         1.856 0.064023 .  
## EggScore2                                                  -0.567 0.571005    
## EggScore3                                                  -1.209 0.227343    
## EggScore4                                                  -2.654 0.008175 ** 
## Test_environmentBlackberry:Original_environmentCherry       0.051 0.959323    
## Test_environmentCherry:Original_environmentCherry          -0.756 0.450084    
## Test_environmentStrawberry:Original_environmentCherry          NA       NA    
## Test_environmentBlackberry:Original_environmentStrawberry   0.454 0.650315    
## Test_environmentCherry:Original_environmentStrawberry          NA       NA    
## Test_environmentStrawberry:Original_environmentStrawberry      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1483 on 568 degrees of freedom
## Multiple R-squared:  0.505,  Adjusted R-squared:  0.4762 
## F-statistic: 17.56 on 33 and 568 DF,  p-value: < 2.2e-16
## F test for SA
(Fratio_int = (anova(m3)[3,2]/anova(m3)[5,2])/(1/anova(m3)[5,1]))
## [1] 64.43531
(pvalue_int = 1 - pf(Fratio_int, 1, anova(m3)[5,1])) #the correct test (see equation D7 in Appendix D of the paper)
## [1] 0.004036806
## Compute Rsquare
(rsq2 <- 1-anova(m2)[5, 3]/((anova(m2)[3, 2]+anova(m2)[5, 2])/(anova(m2)[3, 1]+anova(m2)[5, 1])))
## [1] 0.9642162
## Compare with 5 egg scores
m4 <- lm(asin(sqrt(Emergence_Rate)) ~ Population + Test_environment + SA + Original_environment:Test_environment+EggScoreFive, data=data[data$Generation=="G2",])

## Compare with EggScoreSmall
m5 <- lm(asin(sqrt(Emergence_Rate)) ~ Population + Test_environment + SA + Original_environment:Test_environment+EggScoreSmall, data=data[data$Generation=="G2",])


summary(m4)
## 
## Call:
## lm(formula = asin(sqrt(Emergence_Rate)) ~ Population + Test_environment + 
##     SA + Original_environment:Test_environment + EggScoreFive, 
##     data = data[data$Generation == "G2", ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.42998 -0.09038 -0.00473  0.07960  0.93896 
## 
## Coefficients: (3 not defined because of singularities)
##                                                            Estimate Std. Error
## (Intercept)                                                0.601300   0.055289
## PopulationBlackberry32                                    -0.377346   0.037959
## PopulationBlackberry33                                    -0.055881   0.036355
## PopulationBlackberry34                                    -0.228972   0.053250
## PopulationBlackberry35                                    -0.209601   0.036206
## PopulationBlackberry36                                    -0.107051   0.046078
## PopulationBlackberry37                                    -0.119387   0.038724
## PopulationBlackberry38                                    -0.320919   0.072830
## PopulationBlackberry39                                    -0.107259   0.041814
## PopulationBlackberry40                                     0.012021   0.039828
## PopulationBlackberry43                                    -0.022641   0.055390
## PopulationBlackberry44                                    -0.026557   0.039715
## PopulationBlackberry45                                    -0.236677   0.079817
## PopulationCherry103                                        0.044764   0.049636
## PopulationCherry104                                       -0.240816   0.052751
## PopulationCherry3                                         -0.481162   0.061106
## PopulationCherry47                                        -0.241190   0.042944
## PopulationCherry50                                         0.159812   0.057307
## PopulationCherry51                                         0.278492   0.155026
## PopulationCherry52                                         0.135919   0.083205
## PopulationCherry6                                         -0.360564   0.068615
## PopulationCherry7                                         -0.190717   0.043350
## PopulationStrawberry42                                    -0.007520   0.041778
## PopulationStrawberry44                                    -0.255064   0.044849
## PopulationStrawberry53                                    -0.136571   0.040989
## Test_environmentCherry                                     0.001734   0.047496
## Test_environmentStrawberry                                -0.117162   0.037091
## SA1                                                        0.069558   0.037314
## EggScoreFive2                                             -0.016379   0.028899
## EggScoreFive3                                             -0.036700   0.030340
## EggScoreFive4                                             -0.086037   0.033727
## EggScoreFive5                                             -0.099671   0.047246
## Test_environmentBlackberry:Original_environmentCherry      0.003612   0.048933
## Test_environmentCherry:Original_environmentCherry         -0.043605   0.059561
## Test_environmentStrawberry:Original_environmentCherry            NA         NA
## Test_environmentBlackberry:Original_environmentStrawberry  0.029804   0.064476
## Test_environmentCherry:Original_environmentStrawberry            NA         NA
## Test_environmentStrawberry:Original_environmentStrawberry        NA         NA
##                                                           t value Pr(>|t|)    
## (Intercept)                                                10.876  < 2e-16 ***
## PopulationBlackberry32                                     -9.941  < 2e-16 ***
## PopulationBlackberry33                                     -1.537 0.124831    
## PopulationBlackberry34                                     -4.300 2.01e-05 ***
## PopulationBlackberry35                                     -5.789 1.17e-08 ***
## PopulationBlackberry36                                     -2.323 0.020518 *  
## PopulationBlackberry37                                     -3.083 0.002148 ** 
## PopulationBlackberry38                                     -4.406 1.26e-05 ***
## PopulationBlackberry39                                     -2.565 0.010569 *  
## PopulationBlackberry40                                      0.302 0.762906    
## PopulationBlackberry43                                     -0.409 0.682871    
## PopulationBlackberry44                                     -0.669 0.503964    
## PopulationBlackberry45                                     -2.965 0.003151 ** 
## PopulationCherry103                                         0.902 0.367523    
## PopulationCherry104                                        -4.565 6.13e-06 ***
## PopulationCherry3                                          -7.874 1.75e-14 ***
## PopulationCherry47                                         -5.616 3.06e-08 ***
## PopulationCherry50                                          2.789 0.005470 ** 
## PopulationCherry51                                          1.796 0.072960 .  
## PopulationCherry52                                          1.634 0.102911    
## PopulationCherry6                                          -5.255 2.10e-07 ***
## PopulationCherry7                                          -4.399 1.30e-05 ***
## PopulationStrawberry42                                     -0.180 0.857225    
## PopulationStrawberry44                                     -5.687 2.07e-08 ***
## PopulationStrawberry53                                     -3.332 0.000919 ***
## Test_environmentCherry                                      0.037 0.970893    
## Test_environmentStrawberry                                 -3.159 0.001669 ** 
## SA1                                                         1.864 0.062819 .  
## EggScoreFive2                                              -0.567 0.571108    
## EggScoreFive3                                              -1.210 0.226929    
## EggScoreFive4                                              -2.551 0.011004 *  
## EggScoreFive5                                              -2.110 0.035331 *  
## Test_environmentBlackberry:Original_environmentCherry       0.074 0.941190    
## Test_environmentCherry:Original_environmentCherry          -0.732 0.464411    
## Test_environmentStrawberry:Original_environmentCherry          NA       NA    
## Test_environmentBlackberry:Original_environmentStrawberry   0.462 0.644079    
## Test_environmentCherry:Original_environmentStrawberry          NA       NA    
## Test_environmentStrawberry:Original_environmentStrawberry      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1484 on 567 degrees of freedom
## Multiple R-squared:  0.5051, Adjusted R-squared:  0.4754 
## F-statistic: 17.02 on 34 and 567 DF,  p-value: < 2.2e-16
model.sel(m1, m2, m3, m4, m5)
## Model selection table 
##     (Int) Ppl SA Tst_env Org_env:Tst_env log(Nb_egg+1) EgS ESF ESS df  logLik
## m2 0.9142   +  +       +               +      -0.07308             33 318.394
## m3 0.6020   +  +       +               +                 +         35 312.255
## m5 0.5991   +  +       +               +                         + 37 314.175
## m4 0.6013   +  +       +               +                     +     36 312.318
## m1 0.5794   +  +       +               +                           32 305.067
##      AICc delta weight
## m2 -566.8  0.00      1
## m3 -550.1 16.78      0
## m5 -549.4 17.47      0
## m4 -547.9 18.92      0
## m1 -542.4 24.42      0
## Models ranked by AICc(x)
## Cl= model m5 provides a better description of the data than model m4 (variation in the number of eggs between 50 and 100 is more important than between 150 and 250)

2.5.6 Power Analysis

m2 <- lm(asin(sqrt(Emergence_Rate)) ~ -1 + Population + Test_environment + Original_environment:Test_environment, data=data[data$Generation=="G2",])

simdata <- function(seed=1, nrep = 2, npop = 3, varres=0){
  #print(seed)
  set.seed(seed)
  ## Sample populations
  PopBlackberry <- sort(sample(levels(data$Population)[grep("Blackberry", levels(data$Population))], size=npop, replace = TRUE))
  PopCherry <- sort(sample(levels(data$Population)[grep("Cherry", levels(data$Population))], size=npop, replace = TRUE))
  PopStrawberry <- sort(sample(levels(data$Population)[grep("Strawberry", levels(data$Population))], size=npop, replace = TRUE))
  
  ## Create dataset
  newdata <- expand.grid(Pop_new_name=1:(3*npop), Test_environment=levels(data$Test_environment), Rep=1:nrep)
  newdata$Original_environment <- rep(c("Blackberry", "Cherry", "Strawberry"), each=npop)
  
  ## Add original names to get the right coef in the model
  newdata$Pop_new_name <- as.factor(newdata$Pop_new_name)
  newdata$Population <- newdata$Pop_new_name
  levels(newdata$Population) <- c(PopBlackberry, PopCherry, PopStrawberry)
  
  ## Get design matrix
  mat <- model.matrix(~-1 + Population + Test_environment + Original_environment:Test_environment, data=newdata)
  
  ## Check that the order of the coefficients is the same
  data.frame(colnames(mat), names(coef(m2))[names(coef(m2))%in%colnames(mat)])
  
  coefsim <- coef(m2)[names(coef(m2))%in%colnames(mat)]
  coefsim <- ifelse(is.na(coefsim), 0, coefsim)
  newdata$resp <- mat%*%coefsim + rnorm(n=nrow(newdata), mean = 0, sd = varres)
  
  m3 <- lm(resp ~ -1 + Population + Test_environment + Original_environment:Test_environment, data=newdata)

  data.frame(coefsim, coef(m3))
  
  newdata$SA <- as.factor(ifelse(newdata$Original_environment==newdata$Test_environment, 1, 0))
  
  m4 <- lm(resp ~ Pop_new_name + Test_environment + SA + Original_environment:Test_environment, data=newdata)
  summary(m4)
  anova(m4)
    
  ## F test for SA
  (Fratio_int = (anova(m4)[3,2]/anova(m4)[4,2])/(1/anova(m4)[4,1]))
  (pvalue_int = 1 - pf(Fratio_int, 1, anova(m4)[4,1])) #the correct test (see equation D7 in Appendix D of the paper)
  return(pvalue_int)
}

## With 6 populations per fruit
simdata(seed=2, nrep = 30, npop = 6)
simpval <- sapply(1:500, simdata, nrep = 10, npop = 6, varres=sigma(m2))
mean(simpval)
hist(simpval)
## Compute power
mean(simpval<0.05)

simpval10rep <- sapply(1:500, simdata, nrep = 10, npop = 3, varres=sigma(m2))
simpval20rep <- sapply(1:500, simdata, nrep = 20, npop = 3, varres=sigma(m2))
simpval25rep <- sapply(1:500, simdata, nrep = 25, npop = 3, varres=sigma(m2))
simpval30rep <- sapply(1:500, simdata, nrep = 30, npop = 3, varres=sigma(m2))
simpval35rep <- sapply(1:500, simdata, nrep = 35, npop = 3, varres=sigma(m2))
simpval40rep <- sapply(1:500, simdata, nrep = 40, npop = 3, varres=sigma(m2))
simpval60rep <- sapply(1:500, simdata, nrep = 60, npop = 3, varres=sigma(m2))
simpval80rep <- sapply(1:500, simdata, nrep = 80, npop = 3, varres=sigma(m2))
simpval100rep <- sapply(1:500, simdata, nrep = 100, npop = 3, varres=sigma(m2))

## Compute power
val <- data.frame(simpval10rep, simpval20rep, simpval25rep, simpval30rep, simpval35rep, simpval40rep, simpval60rep, simpval80rep, simpval100rep)

computepower <- function(x){
  return(mean(x<0.05))
}

statpower <- apply(val, 2, computepower)

plot(c(10, 20, 25, 30, 35, 40, 60, 80, 100), statpower, las=1, xlab="Number of replicates", ylab="Statistical Power")

3 Preference

3.1 Load data

setwd("/Users/rodenico/Documents/Pro/Articles/2021_NGS_EvolExp/Simulations")

data <- read.table("DATACOMPLET_PREF_POPNAT2018.csv", header=TRUE, sep=";")
head(data)
##   Generation Experiment Box     Date_P Original_environment Pop_code  Pop_name
## 1         G0 Plasticity 860 19/09/2018           Blackberry       33 Cimetiere
## 2         G0 Plasticity 860 19/09/2018           Blackberry       33 Cimetiere
## 3         G0 Plasticity 860 19/09/2018           Blackberry       33 Cimetiere
## 4         G0 Plasticity 860 19/09/2018           Blackberry       33 Cimetiere
## 5         G0 Plasticity 860 19/09/2018           Blackberry       33 Cimetiere
## 6         G0 Plasticity 860 19/09/2018           Blackberry       33 Cimetiere
##   Pop_original_name   Population Row Column Test_environment Nb_eggs   Date_C_O
## 1      Blackberry33 Blackberry33   1      1        Cranberry       0 13/12/2018
## 2      Blackberry33 Blackberry33   1      2              Fig       0 13/12/2018
## 3      Blackberry33 Blackberry33   1      3        Raspberry       0 13/12/2018
## 4      Blackberry33 Blackberry33   1      4         Rosehips       0 13/12/2018
## 5      Blackberry33 Blackberry33   2      1             Kiwi       0 13/12/2018
## 6      Blackberry33 Blackberry33   2      2       Strawberry       1 13/12/2018
##   Obs_O
## 1    CD
## 2    CD
## 3    CD
## 4    CD
## 5    CD
## 6    CD
data <- data[data$Generation=="G2",]

#Remove GF
data <- data[data$Test_environment!="GF",]

data$Test_environment <- factor(data$Test_environment)

## Remove WT3
data <- data[data$Population!="WT3",]
data$Original_environment <- factor(data$Original_environment)

## Remove Cherry52 which has a very high mean
#data <- data[data$Population!="Cherry52",]
data$Population <- factor(data$Population)
data$Nb_eggs <- as.numeric(as.character(data$Nb_eggs))

data$Row_Col <- as.factor(paste(data$Row, data$Column, sep="_"))

head(data)
##      Generation Experiment Box     Date_P Original_environment Pop_code
## 2665         G2 Adaptation 999 14/11/2018           Blackberry       45
## 2666         G2 Adaptation 999 14/11/2018           Blackberry       45
## 2667         G2 Adaptation 999 14/11/2018           Blackberry       45
## 2668         G2 Adaptation 999 14/11/2018           Blackberry       45
## 2669         G2 Adaptation 999 14/11/2018           Blackberry       45
## 2670         G2 Adaptation 999 14/11/2018           Blackberry       45
##      Pop_name Pop_original_name   Population Row Column Test_environment
## 2665     Uzes      Blackberry45 Blackberry45   1      1            Grape
## 2666     Uzes      Blackberry45 Blackberry45   1      2       Blackberry
## 2667     Uzes      Blackberry45 Blackberry45   1      3           Cherry
## 2668     Uzes      Blackberry45 Blackberry45   1      4       Strawberry
## 2669     Uzes      Blackberry45 Blackberry45   2      1              Fig
## 2670     Uzes      Blackberry45 Blackberry45   2      2         Rosehips
##      Nb_eggs   Date_C_O Obs_O Row_Col
## 2665       0 29/11/2018    LO     1_1
## 2666       0 29/11/2018    LO     1_2
## 2667      16 29/11/2018    LO     1_3
## 2668      12 29/11/2018    LO     1_4
## 2669       0 29/11/2018    LO     2_1
## 2670       0 29/11/2018    LO     2_2
sort(unique(data$Nb_eggs))
##   [1]   0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17
##  [19]  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35
##  [37]  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53
##  [55]  54  55  56  57  58  59  60  61  62  63  64  65  66  67  68  70  71  72
##  [73]  74  75  76  77  78  79  80  81  82  83  84  85  87  89  90  95  98 100
##  [91] 101 107 110 111 115 124 125 128 130 139 141 159

3.2 3 fruits

3.2.1 Explore data with blue

m0 <- lm(Nb_eggs~Original_environment*Test_environment, data=data[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry",])

summary(m0)
## 
## Call:
## lm(formula = Nb_eggs ~ Original_environment * Test_environment, 
##     data = data[data$Test_environment == "Blackberry" | data$Test_environment == 
##         "Cherry" | data$Test_environment == "Strawberry", ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -30.385 -13.622  -5.925   8.260 134.075 
## 
## Coefficients:
##                                                           Estimate Std. Error
## (Intercept)                                                 24.622      2.196
## Original_environmentCherry                                  -8.526      3.730
## Original_environmentStrawberry                              -2.697      4.079
## Test_environmentCherry                                      -1.204      3.106
## Test_environmentStrawberry                                 -17.684      3.106
## Original_environmentCherry:Test_environmentCherry           15.493      5.275
## Original_environmentStrawberry:Test_environmentCherry        4.204      5.769
## Original_environmentCherry:Test_environmentStrawberry       12.261      5.275
## Original_environmentStrawberry:Test_environmentStrawberry    8.534      5.769
##                                                           t value Pr(>|t|)    
## (Intercept)                                                11.212  < 2e-16 ***
## Original_environmentCherry                                 -2.286  0.02263 *  
## Original_environmentStrawberry                             -0.661  0.50871    
## Test_environmentCherry                                     -0.388  0.69840    
## Test_environmentStrawberry                                 -5.694 2.01e-08 ***
## Original_environmentCherry:Test_environmentCherry           2.937  0.00345 ** 
## Original_environmentStrawberry:Test_environmentCherry       0.729  0.46645    
## Original_environmentCherry:Test_environmentStrawberry       2.324  0.02047 *  
## Original_environmentStrawberry:Test_environmentStrawberry   1.479  0.13963    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.74 on 561 degrees of freedom
## Multiple R-squared:  0.1117, Adjusted R-squared:  0.09901 
## F-statistic: 8.816 on 8 and 561 DF,  p-value: 2.139e-11
## Compute proportions
sumeggsperbox <- tapply(data$Nb_eggs, data$Box, sum)
data$sumeggsperbox <- sumeggsperbox[match(as.character(data$Box), names(sumeggsperbox))]
data$prop <-  data$Nb_eggs/data$sumeggsperbox

## Check for local adaptation pattern in proportion
tapply(data$prop[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry"], list(data$Original_environment[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry"], data$Test_environment[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry"]), mean)
##            Apricot Blackberry Blackcurrant    Cherry Cranberry Fig Grape Kiwi
## Blackberry      NA 0.12654832           NA 0.1309270        NA  NA    NA   NA
## Cherry          NA 0.08179375           NA 0.1952284        NA  NA    NA   NA
## Strawberry      NA 0.12029602           NA 0.1280855        NA  NA    NA   NA
##            Raspberry Rosehips Strawberry Tomato
## Blackberry        NA       NA 0.03872770     NA
## Cherry            NA       NA 0.06151998     NA
## Strawberry        NA       NA 0.08369925     NA
## Check for local adaptation pattern in te nulber of eggs
## By pop
tapply(data$Nb_eggs[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry"], list(data$Population[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry"], data$Test_environment[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry"]), mean)
##              Apricot Blackberry Blackcurrant    Cherry Cranberry Fig Grape Kiwi
## Blackberry31      NA  31.250000           NA 24.250000        NA  NA    NA   NA
## Blackberry32      NA  12.000000           NA 13.307692        NA  NA    NA   NA
## Blackberry33      NA  28.000000           NA 32.583333        NA  NA    NA   NA
## Blackberry34      NA  10.000000           NA  6.333333        NA  NA    NA   NA
## Blackberry35      NA  34.857143           NA 23.571429        NA  NA    NA   NA
## Blackberry36      NA  37.800000           NA 35.200000        NA  NA    NA   NA
## Blackberry37      NA  25.181818           NA 22.818182        NA  NA    NA   NA
## Blackberry38      NA  55.000000           NA 37.000000        NA  NA    NA   NA
## Blackberry39      NA  12.000000           NA 29.714286        NA  NA    NA   NA
## Blackberry40      NA  15.714286           NA  9.000000        NA  NA    NA   NA
## Blackberry43      NA   9.333333           NA 12.000000        NA  NA    NA   NA
## Blackberry44      NA  31.538462           NA 30.846154        NA  NA    NA   NA
## Blackberry45      NA   0.000000           NA 16.000000        NA  NA    NA   NA
## Cherry103         NA  11.000000           NA 23.142857        NA  NA    NA   NA
## Cherry104         NA   8.800000           NA 17.000000        NA  NA    NA   NA
## Cherry3           NA   8.333333           NA 25.000000        NA  NA    NA   NA
## Cherry47          NA  16.076923           NA 33.307692        NA  NA    NA   NA
## Cherry50          NA  35.500000           NA 37.250000        NA  NA    NA   NA
## Cherry52          NA   0.500000           NA 28.500000        NA  NA    NA   NA
## Cherry6           NA  53.250000           NA 57.250000        NA  NA    NA   NA
## Cherry7           NA   9.000000           NA 27.857143        NA  NA    NA   NA
## Strawberry42      NA  21.125000           NA 54.250000        NA  NA    NA   NA
## Strawberry44      NA  16.444444           NA 25.333333        NA  NA    NA   NA
## Strawberry53      NA  24.347826           NA 14.565217        NA  NA    NA   NA
##              Raspberry Rosehips Strawberry Tomato
## Blackberry31        NA       NA 12.2500000     NA
## Blackberry32        NA       NA  2.0000000     NA
## Blackberry33        NA       NA  8.1666667     NA
## Blackberry34        NA       NA  5.3333333     NA
## Blackberry35        NA       NA  4.3571429     NA
## Blackberry36        NA       NA 15.2000000     NA
## Blackberry37        NA       NA  5.0909091     NA
## Blackberry38        NA       NA  5.0000000     NA
## Blackberry39        NA       NA  4.4285714     NA
## Blackberry40        NA       NA  9.5714286     NA
## Blackberry43        NA       NA  5.0000000     NA
## Blackberry44        NA       NA  9.1538462     NA
## Blackberry45        NA       NA 12.0000000     NA
## Cherry103           NA       NA  0.4285714     NA
## Cherry104           NA       NA  3.6000000     NA
## Cherry3             NA       NA 28.6666667     NA
## Cherry47            NA       NA  8.1538462     NA
## Cherry50            NA       NA 17.5000000     NA
## Cherry52            NA       NA  1.5000000     NA
## Cherry6             NA       NA 40.2500000     NA
## Cherry7             NA       NA  7.7142857     NA
## Strawberry42        NA       NA 20.8750000     NA
## Strawberry44        NA       NA  8.4444444     NA
## Strawberry53        NA       NA 11.6521739     NA
 ## By original environment
tapply(data$Nb_eggs[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry"], list(data$Original_environment[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry"], data$Test_environment[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry"]), mean)
##            Apricot Blackberry Blackcurrant   Cherry Cranberry Fig Grape Kiwi
## Blackberry      NA   24.62245           NA 23.41837        NA  NA    NA   NA
## Cherry          NA   16.09615           NA 30.38462        NA  NA    NA   NA
## Strawberry      NA   21.92500           NA 24.92500        NA  NA    NA   NA
##            Raspberry Rosehips Strawberry Tomato
## Blackberry        NA       NA   6.938776     NA
## Cherry            NA       NA  10.673077     NA
## Strawberry        NA       NA  12.775000     NA
tapply(data$Nb_eggs[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry"], list(data$Original_environment[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry"], data$Test_environment[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry"]), var)
##            Apricot Blackberry Blackcurrant   Cherry Cranberry Fig Grape Kiwi
## Blackberry      NA   657.9694           NA 515.5036        NA  NA    NA   NA
## Cherry          NA   381.8533           NA 623.6531        NA  NA    NA   NA
## Strawberry      NA   722.4814           NA 836.4814        NA  NA    NA   NA
##            Raspberry Rosehips Strawberry Tomato
## Blackberry        NA       NA   85.29518     NA
## Cherry            NA       NA  405.36161     NA
## Strawberry        NA       NA  264.28141     NA
range(data$Nb_eggs)
## [1]   0 159
## Check for the presence of negative correlations after accounting for the heterogeneity between populations (i.e. some lay few or a lot of eggs)
m1 <- lm(log(Nb_eggs+1) ~ Population, data=data[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry",])
data$resid <- NA
data$resid[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry"] <- residuals(m1)

meanbypopbytestenv <- as.data.frame(tapply(data$resid[data$Generation=="G2"], list(data$Population[data$Generation=="G2"], data$Test_environment[data$Generation=="G2"]), mean))

## Cherry ~ Blackberry
plot(meanbypopbytestenv$Blackberry, meanbypopbytestenv$Cherry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), xlab="Blackberry", ylab="Cherry", pch=16)
legend("bottomleft", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)

## Strawberry ~ Cherry
plot(meanbypopbytestenv$Cherry, meanbypopbytestenv$Strawberry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), xlab="Cherry", ylab="Strawberry", pch=16)
legend("topright", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)

## Strawberry ~ Blackberry
plot(meanbypopbytestenv$Blackberry~ meanbypopbytestenv$Strawberry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), xlab="Blackberry", ylab="Strawberry", pch=16)
legend("bottomleft", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)

3.2.2 Explore data with BLUPS

m0 <- lmer(log(Nb_eggs+1)~ Population + Test_environment + (0 + Test_environment | Population), data = data[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry",])

## Correlation between fixed effects
plot(coef(m1), fixef(m0))

summary(m0)

BLUP <- ranef(m0)$Population

pairs(BLUP)

## Positive correlation between blups and blues, as expected
plot(BLUP$Test_environmentBlackberry, meanbypopbytestenv$Blackberry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), pch=16)
legend("topright", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)

## Positive correlation betwen blups and blues, as expected
plot(BLUP$Test_environmentCherry, meanbypopbytestenv$Cherry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), pch=16)
legend("topleft", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)

## Negative correlation betwen blups and blues, weird
plot(BLUP$Test_environmentStrawberry, meanbypopbytestenv$Strawberry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), pch=16)
legend("bottomright", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)

3.2.3 Check that model with box random effect is OK for the analyses for analyses

data$SA <- as.factor(ifelse(as.character(data$Original_environment)==as.character(data$Test_environment), 1, 0))
data$Box <- as.factor(data$Box)
data$Obs <- as.factor(1:nrow(data))
npop=3

PopBlackberry <- sort(sample(levels(data$Population)[grep("Blackberry", levels(data$Population))], size=npop, replace = TRUE))
  PopCherry <- sort(sample(levels(data$Population)[grep("Cherry", levels(data$Population))], size=npop, replace = TRUE))
  PopStrawberry <- sort(sample(levels(data$Population)[grep("Strawberry", levels(data$Population))], size=npop, replace = TRUE))
  
m1 <- glmer(Nb_eggs ~ 0 + Original_environment + Test_environment + SA + Original_environment:Test_environment+ (1|Box)+(1|Obs), data=data[(data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry")&data$Population%in%c(PopBlackberry, PopCherry, PopStrawberry),], family="poisson")
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00691162 (tol = 0.002, component 1)
summary(m1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Nb_eggs ~ 0 + Original_environment + Test_environment + SA +  
##     Original_environment:Test_environment + (1 | Box) + (1 |      Obs)
##    Data: 
## data[(data$Test_environment == "Blackberry" | data$Test_environment ==  
##     "Cherry" | data$Test_environment == "Strawberry") & data$Population %in%  
##     c(PopBlackberry, PopCherry, PopStrawberry), ]
## 
##      AIC      BIC   logLik deviance df.resid 
##   1110.4   1143.3   -544.2   1088.4      136 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.11291 -0.38074  0.02317  0.10502  0.46559 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Obs    (Intercept) 1.5308   1.2372  
##  Box    (Intercept) 0.5562   0.7458  
## Number of obs: 147, groups:  Obs, 147; Box, 49
## 
## Fixed effects:
##                                                       Estimate Std. Error
## Original_environmentBlackberry                         1.14738    0.36282
## Original_environmentCherry                             1.65103    0.48325
## Original_environmentStrawberry                         2.34639    0.36725
## Test_environmentCherry                                 0.56461    0.43532
## Test_environmentStrawberry                            -0.73526    0.31725
## SA1                                                    0.52175    0.31750
## Original_environmentCherry:Test_environmentCherry     -0.08221    0.85490
## Original_environmentStrawberry:Test_environmentCherry  0.45414    0.53071
## Original_environmentCherry:Test_environmentStrawberry  0.24455    0.67875
##                                                       z value Pr(>|z|)    
## Original_environmentBlackberry                          3.162 0.001565 ** 
## Original_environmentCherry                              3.417 0.000634 ***
## Original_environmentStrawberry                          6.389 1.67e-10 ***
## Test_environmentCherry                                  1.297 0.194625    
## Test_environmentStrawberry                             -2.318 0.020473 *  
## SA1                                                     1.643 0.100320    
## Original_environmentCherry:Test_environmentCherry      -0.096 0.923393    
## Original_environmentStrawberry:Test_environmentCherry   0.856 0.392154    
## Original_environmentCherry:Test_environmentStrawberry   0.360 0.718623    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             Orgn_B Orgn_C Orgn_S Tst_nC Tst_nS SA1    O_C:T_C O_S:T_
## Orgnl_nvrnC  0.023                                                  
## Orgnl_nvrnS  0.393  0.010                                           
## Tst_nvrnmnC -0.650 -0.008 -0.321                                    
## Tst_nvrnmnS -0.393  0.002 -0.433  0.329                             
## SA1         -0.491 -0.006 -0.437  0.405  0.005                      
## Orgnl_C:T_C  0.504 -0.432  0.322 -0.657 -0.170 -0.576               
## Orgnl_S:T_C  0.262  0.000 -0.261 -0.598  0.030 -0.030  0.316        
## Orgnl_C:T_S  0.180 -0.547  0.201 -0.152 -0.468 -0.001  0.385  -0.014
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00691162 (tol = 0.002, component 1)
m1bis <- glmer(Nb_eggs ~ 0 + Original_environment + Test_environment + (1|Box)+(1|Obs), data=data[(data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry")&data$Population%in%c(PopBlackberry, PopCherry, PopStrawberry),], family="poisson")

summary(m1bis)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Nb_eggs ~ 0 + Original_environment + Test_environment + (1 |  
##     Box) + (1 | Obs)
##    Data: 
## data[(data$Test_environment == "Blackberry" | data$Test_environment ==  
##     "Cherry" | data$Test_environment == "Strawberry") & data$Population %in%  
##     c(PopBlackberry, PopCherry, PopStrawberry), ]
## 
##      AIC      BIC   logLik deviance df.resid 
##   1106.4   1127.4   -546.2   1092.4      140 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.14374 -0.37310  0.02474  0.10761  0.40015 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Obs    (Intercept) 1.6136   1.2703  
##  Box    (Intercept) 0.5291   0.7274  
## Number of obs: 147, groups:  Obs, 147; Box, 49
## 
## Fixed effects:
##                                Estimate Std. Error z value Pr(>|z|)    
## Original_environmentBlackberry   1.2851     0.2959   4.343 1.41e-05 ***
## Original_environmentCherry       1.8466     0.3710   4.978 6.43e-07 ***
## Original_environmentStrawberry   2.6376     0.3051   8.644  < 2e-16 ***
## Test_environmentCherry           0.6224     0.2768   2.248   0.0246 *  
## Test_environmentStrawberry      -0.6923     0.2866  -2.415   0.0157 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             Orgn_B Orgn_C Orgn_S Tst_nC
## Orgnl_nvrnC  0.269                     
## Orgnl_nvrnS  0.305  0.251              
## Tst_nvrnmnC -0.496 -0.406 -0.475       
## Tst_nvrnmnS -0.451 -0.377 -0.462  0.497
simdata <- function(seed=1, nrep = 1000, npop = 3, SA=FALSE){
  set.seed(seed)

  ## Create dataset
  newdata <- expand.grid(Original_environment=levels(data$Original_environment), Test_environment=c("Blackberry", "Cherry", "Strawberry"), Box=1:nrep)
  newdata$SA <- as.factor(ifelse(as.character(newdata$Original_environment)==as.character(newdata$Test_environment), 1, 0))
  newdata$Box <- as.factor(newdata$Box)
  newdata$Obs <- as.factor(1:nrow(newdata))

  if(SA==TRUE){
    ## Get design matrix
    mat <- model.matrix(~0 + Original_environment + Test_environment + SA +  Original_environment:Test_environment , data=newdata)
    m <- m1
    ## Drop Strawberry x Strawberry interaction which is not identifiable
    mat <- mat[, -ncol(mat)]
  }else{
    mat <- model.matrix(~0 + Original_environment + Test_environment , data=newdata)
    m <- m1bis
  }
  
  ## Check that the order of the coefficients is the same
  data.frame(colnames(mat), names(fixef(m))[names(fixef(m))%in%colnames(mat)])
  
  ## Get the coefficient of the right populations
  coefsim <- fixef(m)[names(fixef(m))%in%colnames(mat)]
  coefsim <- ifelse(is.na(coefsim), 0, coefsim)

  ## Add box random effect
  mat2 <- model.matrix(~0 + Box, data=newdata)
  coefrand <- rnorm(n=nrep, sd=sqrt(data.frame(VarCorr(m))$vcov[2]))
  
  ## Simulate data
  newdata$resp <- exp(mat%*%coefsim + mat2%*%coefrand + rnorm(n=nrow(newdata), mean = 0, sd = sqrt(as.data.frame(VarCorr(m))$vcov[2])))
  
  newdata$Nb_eggs <- sapply(newdata$resp, rpois, n=1)
  if(SA==TRUE){
    ## Check fit of interaction model
  m3 <- glmer(Nb_eggs ~ 0 + Original_environment + Test_environment + SA + Original_environment:Test_environment + (1|Box) + (1|Obs), data=newdata, family="poisson")
  
  data.frame(coefsim, fixef(m3))
  }
  
  
  m4 <- lmer(log(Nb_eggs+1) ~ 0 + Original_environment + Test_environment + SA + Original_environment:Test_environment + (1|Box) , data=newdata)
  
  anova(m4)
  ## F test for SA
  (Fratio_mixed = (anova(m4)[3,2]/anova(m4)[4,2])/(1/anova(m4)[4,1]))
  ## Compute P-value
  (pvalue_mixed = 1 - pf(Fratio_mixed, 1, anova(m4)[4,1]))
  
  ## Fit SA model
  m5 <- lm(log(Nb_eggs+1) ~ Original_environment + Test_environment + SA + Original_environment:Test_environment + Box, data=newdata)
  summary(m5)
  anova(m5)
    
  ## F test for SA
  (Fratio_fixed = (anova(m5)[3,2]/anova(m5)[5,2])/(1/anova(m5)[5,1]))
  ## Compute P-value
  (pvalue_fixed = 1 - pf(Fratio_fixed, 1, anova(m5)[5,1]))
  return(c(Fratio_mixed=Fratio_mixed, pvalue_mixed=pvalue_mixed, Fratio_fixed=Fratio_fixed, pvalue_fixed=pvalue_fixed))
}

simpval <- sapply(1:100, simdata, nrep = 30, npop = 3, SA=TRUE)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.209911 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.030907 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.137669 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.181223 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0898904 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00268656 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.24508 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0842592 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0526515 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0970607 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00440259 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0135188 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0193863 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00596128 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.012778 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0459003 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00562554 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0624345 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0265195 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00780001 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0320825 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0383621 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00834037 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0663908 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0281588 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0544856 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0595906 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.048528 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.057708 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0494004 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0388851 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0379331 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0139703 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0831132 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0470658 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0554048 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0607758 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0552343 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0163608 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00481468 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0161518 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00475671 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00931366 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0624666 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0612396 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00323709 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.159664 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00256289 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0713899 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.134106 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00527572 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0245835 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.056269 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.10428 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0534692 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0784539 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0135546 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0571271 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.05802 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0536986 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0403313 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0586634 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.057353 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00429039 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0662628 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0454432 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0245136 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.121999 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0634159 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0576719 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0608196 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00231643 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00679617 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0377365 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0457307 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0580104 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0192054 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0220035 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0328747 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0576076 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0592644 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0615031 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0521504 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0543718 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.238321 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.105513 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.156114 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0157434 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0206585 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0582537 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0117851 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0335173 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0452748 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0288024 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00378826 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0661943 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0520898 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00386292 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0375417 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0219405 (tol = 0.002, component 1)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## Power of 51%
mean(simpval[2,]<0.05)
## [1] 0.12
## Check false positive rate of 5%
simpval <- sapply(1:100, simdata, nrep = 100, npop = 3, SA=FALSE)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## False positive rate of 5%
mean(simpval[2,]<0.05)
## [1] 0.02

3.2.4 Analysis

data$SA <- as.factor(ifelse(as.character(data$Original_environment)==as.character(data$Test_environment), 1, 0))
data$Box <- as.factor(data$Box)

## Model with a box effect
m1 <- lmer(log(Nb_eggs+1) ~ Population + Test_environment + SA + (1|Box) + Original_environment:Test_environment, data=data[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry",])
## fixed-effect model matrix is rank deficient so dropping 3 columns / coefficients
## F test for SA with Box effect
(Fratio_int = (anova(m1)[3,2]/anova(m1)[4,2])/(1/anova(m1)[4,1]))
## [1] 11.54878
(pvalue_int = 1 - pf(Fratio_int, 1, anova(m1)[4,1])) #the correct test (see equation D7 in Appendix D of the paper)
## [1] 0.04251289
## Compute Rsquare
(rsq2 <- 1-anova(m1)[4, 3]/((anova(m1)[3, 2]+anova(m1)[4, 2])/(anova(m1)[3, 1]+anova(m1)[4, 1])))
## [1] 0.7250628
## Model with a box effect and a rwo/col effect
m2 <- lmer(log(Nb_eggs+1) ~ Population + Test_environment + SA + (1|Box) + (1|Row_Col) + Original_environment:Test_environment, data=data[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry",])
## fixed-effect model matrix is rank deficient so dropping 3 columns / coefficients
## F test for SA with Box effect
(Fratio_int = (anova(m2)[3,2]/anova(m2)[4,2])/(1/anova(m2)[4,1]))
## [1] 12.46661
(pvalue_int = 1 - pf(Fratio_int, 1, anova(m2)[4,1])) #the correct test (see equation D7 in Appendix D of the paper)
## [1] 0.03861686
## Compute Rsquare
(rsq2 <- 1-anova(m2)[4, 3]/((anova(m2)[3, 2]+anova(m2)[4, 2])/(anova(m2)[3, 1]+anova(m2)[4, 1])))
## [1] 0.7413783

3.2.5 Power Analysis

m2 <- lmer(log(Nb_eggs+1) ~ -1 + Population + Test_environment + Original_environment:Test_environment + (1|Box) + (1|Row_Col), data=data[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry",])
fixef(m2)

simdata <- function(seed=1, nrep = 2, npop = 3){
  #print(seed)
  set.seed(seed)
  ## Sample populations
  PopBlackberry <- sort(sample(levels(data$Population)[grep("Blackberry", levels(data$Population))], size=npop, replace = TRUE))
  PopCherry <- sort(sample(levels(data$Population)[grep("Cherry", levels(data$Population))], size=npop, replace = TRUE))
  PopStrawberry <- sort(sample(levels(data$Population)[grep("Strawberry", levels(data$Population))], size=npop, replace = TRUE))
  
  ## Create dataset
  newdata <- expand.grid(Box=factor(1:nrep), Pop_new_name=1:(3*npop), Test_environment=rep(c("Blackberry", "Cherry", "Strawberry"), each=4))
  newdata$Pop_new_name <- as.factor(newdata$Pop_new_name)
  newdata$Original_environment <- newdata$Pop_new_name
  
  levels(newdata$Original_environment) <- rep(c("Blackberry", "Cherry", "Strawberry"), each=npop)
  newdata$Box <- as.factor(paste(newdata$Box, newdata$Pop_new_name, sep="_"))
  
  #newdata[newdata$Box==2&newdata$Pop_new_name==1,]
  #newdata[newdata$Box=="1_1",]
  #newdata <- newdata[order(newdata$Box),]
  #newdata[newdata$Original_environment=="Blackberry",]

  nlevels(newdata$Box)
  
  ## Add original names to get the right coef in the model
  newdata$Population <- newdata$Pop_new_name
  levels(newdata$Population) <- c(PopBlackberry, PopCherry, PopStrawberry)
  
  unique(newdata$Population)
  
  newdata$Row_Col <- NA
  for (i in levels(newdata$Box)){
    newdata$Row_Col[newdata$Box==i] <- sample(x=1:12, size=12, replace = FALSE)
  }
  newdata$Row_Col <- as.factor(newdata$Row_Col)
  
  ## Get design matrix
  mat <- model.matrix(~-1 + Population + Test_environment + Original_environment:Test_environment, data=newdata)
  ## Remove levels that are not identifiable
  mat <- mat[, !colnames(mat)%in%c("Test_environmentStrawberry:Original_environmentCherry",     "Test_environmentStrawberry:Original_environmentStrawberry")]

  ## Check that the order of the coefficients is the same
  data.frame(colnames(mat), names(fixef(m2))[names(fixef(m2))%in%colnames(mat)])

  ## Get the coefficient of the right populations
  coefsim <- fixef(m2)[names(fixef(m2))%in%colnames(mat)]
  #coefsim <- ifelse(is.na(coefsim), 0, coefsim)

  ## Add box random effect
  mat2 <- model.matrix(~0 + Box, data=newdata)
  coefrand <- rnorm(n=nrep*npop*3, sd=sqrt(data.frame(VarCorr(m2))$vcov[1]))
  
  ## Add row/col random effect
  mat3 <- model.matrix(~0 + Row_Col, data=newdata)
  coefrand2 <- rnorm(n=12, sd=sqrt(data.frame(VarCorr(m2))$vcov[2]))

  ## Simulate data
  newdata$resp <- mat%*%coefsim + mat2%*%coefrand + mat3%*%coefrand2 + rnorm(n=nrow(newdata), mean = 0, sd = sqrt(as.data.frame(VarCorr(m2))$vcov[2]))

  ## Fit interaction model
  m3 <- lmer(resp ~ -1 + Population + Test_environment + Original_environment:Test_environment + (1|Box) + (1|Row_Col), data=newdata)
  data.frame(coefsim, fixef(m3))
  
  ## Fit SA model assuming 9 different populations
  newdata$SA <- as.factor(ifelse(newdata$Original_environment==newdata$Test_environment, 1, 0))
  m4 <- lmer(resp ~ Pop_new_name + Test_environment + SA + Original_environment:Test_environment + (1|Box) + (1|Row_Col), data=newdata)
  summary(m4)
  anova(m4)
    
  ## F test for SA
  (Fratio_int = (anova(m4)[3,2]/anova(m4)[4,2])/(1/anova(m4)[4,1]))
  ## Compute P-value
  (pvalue_int = 1 - pf(Fratio_int, 1, anova(m4)[4,1]))
  return(pvalue_int)
}

simdata(seed=2, nrep = 30, npop = 3)
simpval <- sapply(1:500, simdata, nrep = 10, npop = 3)
mean(simpval)
## Compute power
mean(simpval<0.05)
hist(simpval)

nsim <- 500
simpval10rep <- sapply(1:nsim, simdata, nrep = 10, npop = 3)
simpval15rep <- sapply(1:nsim, simdata, nrep = 15, npop = 3)
simpval20rep <- sapply(1:nsim, simdata, nrep = 20, npop = 3)
simpval25rep <- sapply(1:nsim, simdata, nrep = 25, npop = 3)
simpval30rep <- sapply(1:nsim, simdata, nrep = 30, npop = 3)
simpval35rep <- sapply(1:nsim, simdata, nrep = 35, npop = 3)
simpval40rep <- sapply(1:nsim, simdata, nrep = 40, npop = 3)
simpval60rep <- sapply(1:nsim, simdata, nrep = 60, npop = 3)
simpval80rep <- sapply(1:nsim, simdata, nrep = 80, npop = 3)
simpval100rep <- sapply(1:nsim, simdata, nrep = 100, npop = 3)

## Compute power
val <- data.frame(simpval10rep, simpval15rep, simpval20rep, simpval25rep, simpval30rep, simpval35rep, simpval40rep, simpval60rep, simpval80rep, simpval100rep)

computepower <- function(x){
  return(mean(x<0.05))
}

statpower <- apply(val, 2, computepower)

plot(c(10, 15, 20, 25, 30, 35, 40, 60, 80, 100), statpower, las=1, xlab="Number of replicates", ylab="Statistical Power", main="Preference (3 fruits)")

3.3 12 fruits

3.3.1 Analysis

data$SA <- as.factor(ifelse(as.character(data$Original_environment)==as.character(data$Test_environment), 1, 0))

m1 <- lmer(log(Nb_eggs+1) ~ Population + Test_environment + SA + Original_environment:Test_environment + (1|Box) + (1|Row_Col), data=data[data$Generation=="G2",])
## fixed-effect model matrix is rank deficient so dropping 3 columns / coefficients
summary(m1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## log(Nb_eggs + 1) ~ Population + Test_environment + SA + Original_environment:Test_environment +  
##     (1 | Box) + (1 | Row_Col)
##    Data: data[data$Generation == "G2", ]
## 
## REML criterion at convergence: 7194.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.74710 -0.70531  0.00879  0.68438  2.89632 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Box      (Intercept) 0.38239  0.6184  
##  Row_Col  (Intercept) 0.03124  0.1768  
##  Residual             1.18303  1.0877  
## Number of obs: 2280, groups:  Box, 190; Row_Col, 12
## 
## Fixed effects:
##                                                              Estimate
## (Intercept)                                                  1.506391
## PopulationBlackberry32                                      -0.702176
## PopulationBlackberry33                                       0.391290
## PopulationBlackberry34                                      -0.538466
## PopulationBlackberry35                                      -0.065216
## PopulationBlackberry36                                       0.727676
## PopulationBlackberry37                                      -0.172517
## PopulationBlackberry38                                       0.188417
## PopulationBlackberry39                                      -0.152217
## PopulationBlackberry40                                      -0.407543
## PopulationBlackberry43                                       0.143136
## PopulationBlackberry44                                       0.415712
## PopulationBlackberry45                                      -1.315642
## PopulationCherry103                                         -0.834328
## PopulationCherry104                                         -0.820710
## PopulationCherry3                                            0.213004
## PopulationCherry47                                          -0.286906
## PopulationCherry50                                          -0.066922
## PopulationCherry52                                          -0.275113
## PopulationCherry6                                            0.951440
## PopulationCherry7                                           -0.364385
## PopulationStrawberry42                                       0.841909
## PopulationStrawberry44                                      -0.164577
## PopulationStrawberry53                                      -0.345964
## Test_environmentBlackberry                                   0.525783
## Test_environmentBlackcurrant                                 0.761040
## Test_environmentCherry                                       1.237440
## Test_environmentCranberry                                   -0.065099
## Test_environmentFig                                          0.479561
## Test_environmentGrape                                        1.149459
## Test_environmentKiwi                                         0.726754
## Test_environmentRaspberry                                    0.671069
## Test_environmentRosehips                                     0.356897
## Test_environmentStrawberry                                   0.013150
## Test_environmentTomato                                       0.330173
## SA1                                                          0.560535
## Test_environmentApricot:Original_environmentCherry           0.558490
## Test_environmentBlackberry:Original_environmentCherry        0.235338
## Test_environmentBlackcurrant:Original_environmentCherry      0.212224
## Test_environmentCherry:Original_environmentCherry           -0.007308
## Test_environmentCranberry:Original_environmentCherry         0.415566
## Test_environmentFig:Original_environmentCherry               0.307003
## Test_environmentGrape:Original_environmentCherry            -0.658191
## Test_environmentKiwi:Original_environmentCherry             -0.146252
## Test_environmentRaspberry:Original_environmentCherry         0.130497
## Test_environmentRosehips:Original_environmentCherry          0.266882
## Test_environmentStrawberry:Original_environmentCherry        0.309021
## Test_environmentApricot:Original_environmentStrawberry       0.418854
## Test_environmentBlackberry:Original_environmentStrawberry    0.474579
## Test_environmentBlackcurrant:Original_environmentStrawberry  0.184944
## Test_environmentCherry:Original_environmentStrawberry       -0.123624
## Test_environmentCranberry:Original_environmentStrawberry    -0.006736
## Test_environmentFig:Original_environmentStrawberry          -0.121108
## Test_environmentGrape:Original_environmentStrawberry        -0.685244
## Test_environmentKiwi:Original_environmentStrawberry         -0.409178
## Test_environmentRaspberry:Original_environmentStrawberry     0.156703
## Test_environmentRosehips:Original_environmentStrawberry      0.256809
##                                                             Std. Error t value
## (Intercept)                                                   0.271768   5.543
## PopulationBlackberry32                                        0.311642  -2.253
## PopulationBlackberry33                                        0.316550   1.236
## PopulationBlackberry34                                        0.469519  -1.147
## PopulationBlackberry35                                        0.307372  -0.212
## PopulationBlackberry36                                        0.395370   1.840
## PopulationBlackberry37                                        0.322254  -0.535
## PopulationBlackberry38                                        0.735595   0.256
## PopulationBlackberry39                                        0.358934  -0.424
## PopulationBlackberry40                                        0.358934  -1.135
## PopulationBlackberry43                                        0.469519   0.305
## PopulationBlackberry44                                        0.311642   1.334
## PopulationBlackberry45                                        0.735595  -1.789
## PopulationCherry103                                           0.401009  -2.081
## PopulationCherry104                                           0.433927  -1.891
## PopulationCherry3                                             0.502417   0.424
## PopulationCherry47                                            0.359298  -0.799
## PopulationCherry50                                            0.460805  -0.145
## PopulationCherry52                                            0.576702  -0.477
## PopulationCherry6                                             0.460805   2.065
## PopulationCherry7                                             0.355601  -1.025
## PopulationStrawberry42                                        0.398242   2.114
## PopulationStrawberry44                                        0.389764  -0.422
## PopulationStrawberry53                                        0.345523  -1.001
## Test_environmentBlackberry                                    0.328650   1.600
## Test_environmentBlackcurrant                                  0.155951   4.880
## Test_environmentCherry                                        0.155700   7.948
## Test_environmentCranberry                                     0.155495  -0.419
## Test_environmentFig                                           0.155674   3.081
## Test_environmentGrape                                         0.155742   7.381
## Test_environmentKiwi                                          0.155801   4.665
## Test_environmentRaspberry                                     0.155733   4.309
## Test_environmentRosehips                                      0.156020   2.288
## Test_environmentStrawberry                                    0.156214   0.084
## Test_environmentTomato                                        0.155728   2.120
## SA1                                                           0.289287   1.938
## Test_environmentApricot:Original_environmentCherry            0.264159   2.114
## Test_environmentBlackberry:Original_environmentCherry         0.421525   0.558
## Test_environmentBlackcurrant:Original_environmentCherry       0.264235   0.803
## Test_environmentCherry:Original_environmentCherry             0.359791  -0.020
## Test_environmentCranberry:Original_environmentCherry          0.264373   1.572
## Test_environmentFig:Original_environmentCherry                0.264138   1.162
## Test_environmentGrape:Original_environmentCherry              0.264271  -2.491
## Test_environmentKiwi:Original_environmentCherry               0.264306  -0.553
## Test_environmentRaspberry:Original_environmentCherry          0.264269   0.494
## Test_environmentRosehips:Original_environmentCherry           0.264409   1.009
## Test_environmentStrawberry:Original_environmentCherry         0.264158   1.170
## Test_environmentApricot:Original_environmentStrawberry        0.289341   1.448
## Test_environmentBlackberry:Original_environmentStrawberry     0.501087   0.947
## Test_environmentBlackcurrant:Original_environmentStrawberry   0.289304   0.639
## Test_environmentCherry:Original_environmentStrawberry         0.289450  -0.427
## Test_environmentCranberry:Original_environmentStrawberry      0.289217  -0.023
## Test_environmentFig:Original_environmentStrawberry            0.289228  -0.419
## Test_environmentGrape:Original_environmentStrawberry          0.288957  -2.371
## Test_environmentKiwi:Original_environmentStrawberry           0.289304  -1.414
## Test_environmentRaspberry:Original_environmentStrawberry      0.289345   0.542
## Test_environmentRosehips:Original_environmentStrawberry       0.289024   0.889
## 
## Correlation matrix not shown by default, as p = 57 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 3 columns / coefficients
anova(m1)
## Analysis of Variance Table
##                                       npar  Sum Sq Mean Sq F value
## Population                              23  90.908  3.9525  3.3410
## Test_environment                        11 290.703 26.4276 22.3388
## SA                                       1  18.790 18.7904 15.8832
## Test_environment:Original_environment   21  55.756  2.6551  2.2443
## F test for SA with Box effect
(Fratio_int = (anova(m1)[3,2]/anova(m1)[4,2])/(1/anova(m1)[4,1]))
## [1] 7.077208
(pvalue_int = 1 - pf(Fratio_int, 1, anova(m1)[4,1])) #the correct test (see equation D7 in Appendix D of the paper)
## [1] 0.01464063
## Compute Rsquare
(rsq2 <- 1-anova(m1)[4, 3]/((anova(m1)[3, 2]+anova(m1)[4, 2])/(anova(m1)[3, 1]+anova(m1)[4, 1])))
## [1] 0.2164463

3.3.2 Power Analysis

m2 <- lmer(log(Nb_eggs+1) ~ -1+Population + Test_environment + Original_environment:Test_environment + (1|Box) + (1|Row_Col), data=data[data$Generation=="G2",])

fixef(m2)
hist(log(data$Nb_eggs[data$Generation=="G2"]+1))

simdata <- function(seed=1, nrep = 2, npop = 3){
  set.seed(seed)
  ## Sample populations
  PopBlackberry <- sort(sample(levels(data$Population)[grep("Blackberry", levels(data$Population))], size=npop, replace = TRUE))
  PopCherry <- sort(sample(levels(data$Population)[grep("Cherry", levels(data$Population))], size=npop, replace = TRUE))
  PopStrawberry <- sort(sample(levels(data$Population)[grep("Strawberry", levels(data$Population))], size=npop, replace = TRUE))
  
  ## Create dataset
  newdata <- expand.grid(Pop_new_name=1:(3*npop), Test_environment=levels(data$Test_environment), Box=factor(1:nrep))
  newdata$Original_environment <- rep(c("Blackberry", "Cherry", "Strawberry"), each=npop)
  
  
  ## Add original names to get the right coef in the model
  newdata$Pop_new_name <- as.factor(newdata$Pop_new_name)
  newdata$Population <- newdata$Pop_new_name
  levels(newdata$Population) <- c(PopBlackberry, PopCherry, PopStrawberry)
  newdata$Box <- as.factor(paste(newdata$Box, newdata$Pop_new_name, sep="_"))
  #newdata[newdata$Pop_new_name=="1",]
  #newdata[order(newdata$Box),]
  #nlevels(newdata$Box)
  unique(newdata$Population)
  
  newdata$Row_Col <- NA
  for (i in levels(newdata$Box)){
    newdata$Row_Col[newdata$Box==i] <- sample(x=1:12, size=12, replace = FALSE)
  }
  newdata$Row_Col <- as.factor(newdata$Row_Col)
  
  ## Get design matrix
  mat <- model.matrix(~-1 + Population + Test_environment + Original_environment:Test_environment, data=newdata)
  ## Remove levels that are not identifiable
  mat <- mat[, !colnames(mat)%in%colnames(mat)[!colnames(mat)%in%names(fixef(m2))]]
  
  ## Check that the order of the coefficients is the same
  data.frame(colnames(mat), names(fixef(m2))[names(fixef(m2))%in%colnames(mat)])
  
  coefsim <- fixef(m2)[names(fixef(m2))%in%colnames(mat)]
  #coefsim <- ifelse(is.na(coefsim), 0, coefsim)
  
  ## Add box random effect
  mat2 <- model.matrix(~0 + Box, data=newdata)
  coefrand <- rnorm(n=nrep*npop*3, sd=sqrt(data.frame(VarCorr(m2))$vcov[1]))
  
  ## Add row/col random effect
  mat3 <- model.matrix(~0 + Row_Col, data=newdata)
  coefrand2 <- rnorm(n=12, sd=sqrt(data.frame(VarCorr(m2))$vcov[2]))
  
  ## Simulate data
  newdata$resp <- mat%*%coefsim + mat2%*%coefrand + mat3%*%coefrand2 + rnorm(n=nrow(newdata), mean = 0, sd = sqrt(as.data.frame(VarCorr(m2))$vcov[3]))

  m3 <- lmer(resp ~ -1 + Population + Test_environment + Original_environment:Test_environment + (1|Box) + (1|Row_Col), data=newdata)
  
  data.frame(coefsim, fixef(m3))
  
  newdata$SA <- as.factor(ifelse(newdata$Original_environment==newdata$Test_environment, 1, 0))
  
  m4 <- lmer(resp ~ Pop_new_name + Test_environment + SA + Original_environment:Test_environment + (1|Box) + (1|Row_Col), data=newdata)
  summary(m4)
  anova(m4)
    
  ## F test for SA
  (Fratio_int = (anova(m4)[3,2]/anova(m4)[4,2])/(1/anova(m4)[4,1]))
  (pvalue_int = 1 - pf(Fratio_int, 1, anova(m4)[4,1])) #the correct test (see equation D7 in Appendix D of the paper)
  return(pvalue_int)
}

simdata(seed=2, nrep = 30, npop = 3)
simpval <- sapply(1:500, simdata, nrep = 10, npop = 3)
mean(simpval)
## Compute power
mean(simpval<0.05)
hist(simpval)

nsim <- 500
simpval10rep <- sapply(1:nsim, simdata, nrep = 10, npop = 3)
simpval15rep <- sapply(1:nsim, simdata, nrep = 15, npop = 3)
simpval20rep <- sapply(1:nsim, simdata, nrep = 20, npop = 3)
simpval25rep <- sapply(1:nsim, simdata, nrep = 25, npop = 3)
simpval30rep <- sapply(1:nsim, simdata, nrep = 30, npop = 3)
simpval35rep <- sapply(1:nsim, simdata, nrep = 35, npop = 3)
simpval40rep <- sapply(1:nsim, simdata, nrep = 40, npop = 3)
simpval60rep <- sapply(1:nsim, simdata, nrep = 60, npop = 3)
simpval80rep <- sapply(1:nsim, simdata, nrep = 80, npop = 3)
simpval100rep <- sapply(1:nsim, simdata, nrep = 100, npop = 3)

## Compute power
val <- data.frame(simpval10rep, simpval15rep, simpval20rep, simpval25rep, simpval30rep, simpval35rep, simpval40rep, simpval60rep, simpval80rep, simpval100rep)

computepower <- function(x){
  return(mean(x<0.05))
}

statpower <- apply(val, 2, computepower)

plot(c(10, 15, 20, 25, 30, 35, 40, 60, 80, 100), statpower, las=1, xlab="Number of replicates", ylab="Statistical Power", main="Preference (12 fruits)")